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Abstract 


The  first  completely  physical  electro-thermal  model  is  presented  that  is  capable  of  describing  the 
large  signal  performance  of  MESFET-  and  HEMT-based,  high  power  microwave  and  millimeter  wave 
monolithic  and  hybrid  ICs,  on  timescales  suitable  for  CAD.  The  model  includes  the  effects  of  self¬ 
heating  and  mutual  thermal  interaction  on  active  device  performance  with  full  treatment  of  all  thermal 
non  linearities.  The  electrical  description  is  provided  by  the  rapid  quasi-2D  Leeds  Physical  Model 
and  the  steady-state  global  thermal  description  is  provided  by  a  highly  accurate  and  computationally 
inexpensive  analytical  thermal  resistance  matrix  approach.  The  order  of  the  global  thermal  resistance 
matrix  describing  3-dimensional  heat  flow  in  complex  systems,  is  shown  to  be  determined  purely  by 
the  number  of  active  device  elements,  not  the  level  of  internal  device  structure.  Thermal  updates 
in  the  necessarily  iterative,  fully  coupled  electro-thermal  solution,  therefore  reduce  to  small  matrix 
multiplications  implying  orders  of  magnitude  speed-up  compared  to  the  use  of  full  numerical  thermal 
solutions  capable  of  comparable  levels  of  detail  and  accuracy. 

KEYWORDS:  electrothermal,  thermal,  modelling,  MMIC,  power  transistor. 
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1  Introduction 


The  importance  of  thermal  effects  in  the  description  of  high  power  device  and  circuit  performance  has  long 
been  recognised,  and  much  work  has  been  done  in  modelling  the  effects  of  self-heating  and  mutual  thermal 
interaction  on  device  and  circuit  performance  [1]— [48] .  If  electro-thermal  models  are  to  be  of  value  for 
predictive  electro-thermal  design  studies,  based  purely  on  specified  material  and  structural  information 
rather  than  on  prior  experimental  characterisation,  then  fully  physical  (as  opposed  to  equivalent  circuit  or 
physics-based)  coupled  models  are  required.  However,  although  physical  thermal  simulations  have  been 
performed  for  complex  device  structures,  relatively  little  work  has  been  presented  describing  explicit 
coupling  of  physical  thermal  and  large  signal  electrical  models.  The  reason  for  the  absence  of  such 
coupled  calculations  arises  from  three  main  sources.  Firstly,  fully  physical  device  models  with  large 
signal  electrical  capability  are  still  not  commonplace.  Secondly,  although  the  basic  thermal  issues  are 
well  understood,  thermal  solvers  available  in  other  areas  of  study,  such  as  mechanical  engineering,  are 
not  readily  adaptable  to  electronic  engineering  problems.  Thirdly,  the  coupled  electro-thermal  problem 
is  intrinsically  non  linear  implying  the  need  for  iterative  solutions,  and  existing  physical  thermal  and 
electrical  solvers  are  generally  too  slow  to  generate  an  iterative  coupled  solution  sufficiently  fast  for 
electro- thermal  design  studies. 

Of  the  existing  thermal  models  with  the  potential  for  coupled  electro-thermal  simulations,  that  of 
Bonani  et  al  for  steady-state  simulation  of  heatsink  mounted  MMICs  is  most  comprehensive  [l]-[3]. 
This  model  describes  arbitrary  surface  metallisation,  vias  and  partial  substrate  thinning  for  multi-finger 
MMICs.  Thermal  solutions  based  on  this  hybrid  Green’s  function  finite  element  approach  typically  take 
~1  hour  on  a  medium-sized  workstation  [3].  Webb  has  performed  thermal  simulations  for  both  steady- 
state  and  thermal  transient  cases  [4]— [8] ,  and  typical  steady-state  thermal  simulations  take  ~l/2  hour 
[8]  on  a  PC  using  an  optimised  finite  difference  scheme.  Dorkel  et  al  have  performed  steady-state  and 
transient  thermal  simulations  taking  just  a  few  minutes  but  these  were  for  simplified  multilayer  structures, 
based  on  a  2-port  network  (or  2  x  2  transfer  matrix)  formulation  [9].  All  of  these  calculations  assumed 
adiabatic  surface  boundary  conditions.  Newtonian  boundary  conditions  have  been  imposed  for  use  at 
the  circuit  board  level  where  surface  areas  are  large,  for  instance  by  Ellison  who  performed  analytical 
Green’s  function  calculations  for  simple  multilayer  structures  [10,  11],  Electro-thermal  simulations  have 
been  performed  by  Anholt  for  HBTs,  MESFETs  and  HEMTs  [12]— [15]  and  compared  against  simple 
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analytical  expressions  for  thermal  resistances  [15].  Ghione  and  Naldi  performed  fully  physical,  self- 
consistent  electro-thermal  calculations  for  MESFETs  in  MMICs  [16].  The  boundary  conditions  for  their 
physical  model  were  provided  by  a  thermal  resistance  formulation,  with  the  thermal  resistances  obtained 
approximately  using  a  conformal  mapping  technique. 

It  has  been  stated  repeatedly  [2,  9]  that  the  thermal  resistance  method  is  unsuitable  for  accurate 
description  of  detailed  device  structures,  but  this  assertion  depends  on  a  particular  conception  of  the 
thermal  resistance  approach  as  being  fundamentally  approximate.  This  paper  shows  that  the  thermal 
resistance  matrix  concept  [16]— [19]  can  be  developed  as  an  essentially  analytically  exact  description  of 
solutions  of  the  heat  diffusion  equation  at  selected  points  within  a  complex  3-dimensional  volume.  Tem¬ 
perature  updates  in  the  necessarily  iterative  coupled  electro- thermal  solution,  based  on  this  resistance 
matrix  approach,  reduce  to  small  matrix  multiplications  with  the  order  of  the  matrices  determined  solely 
by  the  number  of  active  device  elements.  Such  solutions  are  therefore  orders  of  magnitude  faster  than 
solutions  based  on  finite  element  or  finite  difference  thermal  descriptions,  whilst  achieving  the  same  levels 
of  detail  and  accuracy. 

The  thermal  resistance  matrix  description  presented  here  is  applicable  from  individual  MMICs  through 
hybrid  MICs,  RFICs  and  MMIC  grid  arrays  to  circuit  board  level,  where  the  number  of  active  device 
elements  is  ~10  to  several  xlO3.  Electro-thermal  simulations  for  VLSI  design  have  been  presented 
previously  [20].  However,  use  of  the  thermal  resistance  matrix  approach  presented  here,  in  an  unmodifed 
form,  with  the  order  of  the  dense  thermal  resistance  matrix  given  by  the  number  of  individual  active 
devices  typical  in  VLSI  or  ULSI,  would  clearly  be  of  no  advantage  compared  to  standard  numerical 
solutions  of  the  heat  diffusion  equation  or  to  sparse  matrix  thermal  network  formulations.  The  thermal 
resistance  matrix  description  given  in  this  paper  is  therefore  not  intended  to  apply  directly  to  the  VLSI 
case. 

The  global  thermal  resistance  matrix  model  described  here  is  given  in  the  context  of  coupled  electro¬ 
thermal  modelling.  However,  it  is  applicable  to  any  thermal  problem  that  can  be  well  described  by  the 
temperatures  in  the  vicinities  of  a  relatively  small  number  of  localised  regions,  rather  than  requiring 
temperature  information  throughout  the  whole  body  of  a  complex  volume.  Use  of  the  Kirchhoff  trans¬ 
formation,  combined  with  use  of  the  resistance  matrix  for  (generally  non  linear)  interface  matching  of 
analytical  subsystem  solutions,  is  more  generally  applicable  to  any  physical  system  described  by  Laplace’s 
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equation 


v.  =  0  (1) 

providing  the  system  can  be  sufficiently  accurately  represented  by  rectangular  subsystems,  with  piecewise 
discretisation  of  ip  and  its  derivative,  ip' ,  at  the  interfaces. 

The  boundary  element  method  replaces  3-dimensional  numerical  finite  element  or  finite  difference  de¬ 
scriptions,  given  by  large  sparse  matrices,  with  a  2-dimensional  surface  description  given  by  smaller  dense 
matrices,  by  using  a  Green’s  function  technique.  In  a  similar  fashion  the  analytical  thermal  resistance 
matrix  approach  reduces  the  thermal  description  to  that  of  very  small  dense  matrices,  describing  only 
limited  2-dimensional  regions  at  surfaces  and  discretised  interfaces  by  direct  solution  of  the  heat  diffusion 
equation.  All  explicit  reference  to  redundant  temperature  information  at  discretised  volume  and  surface 
nodes  is  eliminated.  This  thermal  resistance  matrix  approach  can  equally  be  viewed  as  a  generalised  finite 
element  or  finite  volume  approach,  in  which  the  finite  elements  are  not  primitive  volumes  with  tempera¬ 
tures  between  nodes  described  by  low  order  polynomial  interpolation,  but  complete  thermal  subsystems 
with  internal  temperatures  given  by  full  analytical  solutions  of  the  heat  diffusion  equation.  Construction 
of  a  global  thermal  resistance  matrix  then  consists  of  matrix  manipulations  on  resistance  matrices  given 
by  simple  analytical  expressions  for  thermal  subsystems. 

A  key  feature  of  this  thermal  resistance  matrix  approach  is  that  the  global  thermal  resistance  matrix 
only  has  to  be  evaluated  once,  for  the  thermal  steady-state  case  treated  here,  prior  to  the  coupled 
electro-thermal  solution.  Construction  time  for  the  thermal  resistance  matrices  has  no  impact  on  coupled 
electro-thermal  run  time.  The  global  thermal  resistance  matrix  is  itself  constructed  rapidly  by  simple 
matrix  manipulations  on  analytical  expressions  for  thermal  subelements,  such  as  MMIC  dies  or  sections  of 
surface  metallisation,  each  evaluated  purely  in  terms  of  subsystem  layout  and  material  parameters.  The 
global  resistance  matrix  describes  totally  arbitrary  device  layouts  (based  on  simple  rectangular  structures) 
in  3  dimensions  and  once  constructed  can  be  used  repeatedly  in  coupled  electro-thermal  simulations. 

Models  containing  descriptions  of  self-heating  effects  on  the  performance  of  a  range  of  devices  have 
been  described  previously,  e.g.  [17],  [21]— [23],  The  global  thermal  model  presented  here  is  intended  to  be 
readily  compatible  with  any  thermally  self-consistent  device  model,  without  requiring  core  changes.  It 
should  then  supply  a  fully  physical  and  highly  accurate  description  of  mutual  thermal  interaction  between 
active  devices  at  chip  level  and  above. 


If 


The  paper  provides  an  overview  of  the  analytical  thermal  resistance  matrix  approach,  demonstrating 
how  it  can  achieve  the  same  level  of  accuracy  and  detail  as  finite  element  or  finite  difference  descriptions 
such  as  those  of  Bonani  or  Webb,  whilst  reducing  the  global  thermal  description  to  small  dense  matrix 
multiplications.  This  reduction  implies  rapid  temperature  updates  in  the  coupled  electro-thermal  model 
allowing  fully  physical  solution  on  timescales  suitable  for  CAD.  The  global  thermal  description  combines 
analytical  solutions  of  the  heat  diffusion  equation  for  thermal  subsystems  with  hierarchical  interface 
matching  of  thermal  resistance  matrices  for  complex  3-dimensional  systems.  It  is  shown  that  this  approach 
allows  accurate  treatment  of  all  thermal  non  linearities  arising  from  temperature  dependent  material 
constants  and  surface  fluxes. 

The  paper  is  structured  as  follows.  The  next  section  contains  a  brief  outline  of  the  Leeds  Physical 
Model  (LPM)  which  provides  the  large  signal  active  device  description  required  in  the  coupled  electro¬ 
thermal  model.  This  is  followed  by  a  statement  of  the  coupled  electro- thermal  problem  and  an  indication 
of  the  key  issues  in  the  global  thermal  description.  The  motivation  for  the  thermal  resistance  matrix 
approach  to  the  global  thermal  problem  is  then  indicated.  This  is  followed  by  description  of  the  analytical 
thermal  steady-state  solutions  for  the  MMIC  and  for  other  subsystem  components.  Simple  analytical 
expressions  are  presented  for  the  thermal  resistance  matrices  of  entire  MMICs  bearing  several  multi-finger 
power  transistors.  A  linear  ‘radiation’  boundary  condition  is  imposed  to  describe  convective  and  radiative 
surface  losses  from  large  area  substrates.  A  hierarchical  approach  is  presented  for  the  resistance  matrix 
description  of  multilayered  structures  such  as  MMICs  with  superstructure,  e.g.  surface  metallisation  and 
air  bridges,  or  for  multiple  MMICs  mounted  on  a  common  substrate.  A  double  Fourier  series  finite  element 
technique  is  outlined  for  the  construction  of  the  thermal  resistance  matrices  of  ICs  with  totally  arbitrary 
surface  metallisation  and  full  3-dimensional  heat  flow.  An  original  analytical  solution  is  presented  for  the 
3-dimensional  temperature  distribution  and  corresponding  thermal  resistance  matrix  in  a  hybrid  MIC  or 
MMIC  with  any  number  of  arbitrarily  placed  vias  of  arbitrary  cross-section,  or  with  partial  substrate 
thinning.  Finally,  the  paper  is  summarised  and  conclusions  drawn. 

2  The  Leeds  Physical  Model 

The  Leeds  Physical  Model  [22] ,  [49]— [51]  is  a  quasi-2D  physical  model  of  MESFETs  and  HEMTs  including 
the  effects  of  self-heating.  It  is  a  sophisticated  CAD  tool  and  has  been  released  commercially  as  part  of 
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HP-EEsof’s  MDS  software. 


The  quasi-2D  approximation  is  based  on  the  observation,  from  full  2-dimensional  simulations,  that 
carrier  transport  is  essentially  1-dimensional  and  driven  by  the  component  of  electric  field  along  the  device 
channel.  The  model  solves  a  set  of  hydrodynamic  equations  obtained  from  moments  of  the  Boltzmann 
equation: 


Qy~l  _ _ .  y  v 

—  +  V  ■  (nv) 

dv  _ 

-+„-v„ 

dw  — 

—  + v • Vw 
ot 


ffl 

CLpL~dt 


=  0, 


-t-E  - 


2  V(nw)  +  V(nv 2)  -  — , 


3  m*n 


3  n 

/  2\ 
nv(w  — 2~v  ) 


-Iv.q-Zz*, 

n  rw 


=  V  •  (kl' VT)  +  J  •  E9 


(2) 


in  the  direction  parallel  to  the  heterointerface  (suitably  reduced  for  rapid  solution)  [49].  Here  n  is  carrier 
density,  v  is  velocity,  w  is  energy,  m*  is  effective  mass,  E  is  electric  field,  rp,  tw  are  relaxation  lifetimes, 
T  is  lattice  temperature,  cL  is  heat  capacity,  pl  is  density,  t  is  time,  kl  is  thermal  conductivity  and 
J  is  current  density.  This  energy  transport  model  then  consists  of  continuity,  momentum  and  energy 
conservation,  and  heat  diffusion  equations. 

The  LPM  solves  Poisson  and  Schrodinger  equations  self-consistently  in  the  direction  normal  to  the 
heterointerface,  to  describe  charge  control,  Figure  1.  The  charge  control  information  is  calculated  prior 
to  calculation  of  in-plane  transport  and  is  stored  as  a  look-up  table. 

The  quasi-2D  model  includes  a  full  description  of  the  device  cross-section,  by  describing  charge  con¬ 
servation  in  the  vicinity  of  the  device  channel  via  a  series  of  Gaussian  boxes,  Figure  2. 

The  Leeds  Physical  Model  incorporates  the  effects  of  temperature  on  device  performance  by  use  of  a 
temperature  dependent  low  field  mobility, 

P  =  m(T0)(^)",  (3) 


where  T  is  the  device  channel  temperature  and  n  =  2.3  [52].  Velocity-field  characteristics  are  obtained 
from  Monte  Carlo  simulations  and  stored  as  a  simple  parameterisation. 

The  LPM  provides  a  physical  large  signal  description  of  device  performance,  as  illustrated  by  Figure 
3,  which  shows  calculated  time  and  frequency  domain  response  to  a  sinusoidal  large  signal  input  [53]. 
Output  waveform  distortion  and  generation  of  harmonics  differing  from  the  fundamental  are  clear. 
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The  LPM  forms  the  core  of  the  coupled  electro-thermal  solution  presented  in  this  paper.  The  coupled 
electro- thermal  problem  is  now  described. 

3  The  Coupled  Electro- Thermal  Problem 

The  non  linear  coupled  electro-thermal  problem  can  be  described  by  the  equation 

ei(P1,...,FM)-ei{Pi)  =  0  for  *  =  1, M  (4) 

for  self-consistent  solution  of  active  device  temperatures  6t  and  power  dissipations  Pi,  in  device  channels 
i  =  1 

The  first  term  states  that  the  channel  temperature  of  the  ith  active  device  is  a  (non  linear)  function 
of  the  power  dissipations  Pj  of  all  heat  dissipating  elements,  j  =  as  described  by  the  global 

thermal  model.  The  second  term  states  that  the  ith  channel  temperature  has  a  unique  relation  to  the 
corresponding  power  dissipation  Pi  of  the  ith  active  element  as  determined  by  the  LPM. 

Simultaneous  imposition  of  the  non  linear  description  of  heat  flow  through  the  MMIC  or  hybrid  MIC 
(with  surface  heat  loss),  and  the  non  linear  electrical  description  of  transistor  action,  determine  uniquely 
the  temperatures  and  power  dissipations  in  the  vicinities  of  the  active  device  channels.  These  in  turn, 
determine  full  electrical  solutions,  DC  or  large  signal  RF.  The  above  description  applies  equally  well  to 
both  the  thermal  steady-state  and  to  the  time-dependent  case  in  Laplace  transform  space. 

Details  of  time-dependent  electrical  behaviour  are  decoupled  from  thermal  behaviour  by  the  fact  that 
the  inverse  lattice  thermal  time  constant  is  orders  of  magnitude  different  from  GHz  electrical  frequencies. 
The  large  signal  RF  case  is  therefore  treated  by  time  averaging  the  electrical  signal  (to  describe  effects  such 
as  AC  cooling).  Only  the  steady-state  thermal  case  is  treated  here.  Thermal  resistance  matrix  treatment 
of  thermal  transients  due  to  turn-on,  or  due  to  low  frequency  pulsed  operation,  will  be  described  elsewhere 
[24]- 


3.1  Coupling  the  Leeds  Physical  Model 

The  manner  in  which  the  LPM  enters  the  coupled  solution  is  easily  illustrated  for  the  DC  case.  The  LPM 
is  invoked  as  a  simple  subroutine  call  for  each  active  device  element 

Ids  =  Ids(Vds,  Vgs\  T).  (5) 
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This  device  element  will  typically  be  a  single  source-gate-drain  subunit,  or  a  fraction  of  such  a  unit 
divided  along  the  gate  width. 

Supplying  bias  point  (Vds,  Vgs),  where  Vds  is  drain-source  voltage  and  Vgs  is  gate-source  voltage, 
and  with  active  device  channel  temperature  given  by  T  (assumed  uniform  along  the  channel  for  simplicity 
of  description),  the  LPM  returns  drain-source  current  IDS .  The  active  device  power  dissipation  P  is  then 
given  by  the  approximate  relation, 

P  fit  IdsVds •  (6) 

P  is  the  time  average  value  in  the  time  dependent  case. 

A  single,  uniform  channel  temperature  is  assumed  in  the  current  coupled  electro-thermal  solution, 
but  this  is  not  a  necessary  condition  for  solution.  The  LPM  returns  the  power  dissipation  P  =  E.J  as 
a  function  of  field,  E)  and  current  density,  J,  along  the  device  channel.  By  subdividing  the  gate  along 
its  length  the  resistance  matrix  approach  can  describe  a  varying  channel  temperature,  with  consequent 
improvement  in  transistor  breakdown  modelling  [23]  and  without  the  need  to  artificially  displace  the  gate 
in  order  to  describe  correctly  the  position  of  the  peak  temperature  in  the  vicinity  of  a  gate  finger.  As 
the  LPM  returns  E.J ,  the  extent  of  the  device  heating  is  determined  accurately  by  the  physical  model. 
No  approximation  is  required  to  provide  an  effective  gate  length  for  comparison  against  experiment  [16]. 
The  position  dependent  power  dissipation  in  the  LPM  is  finely  discretised,  but  the  discretisation  used 
for  the  thermal  description,  along  the  source-gate-drain  length,  can  be  much  coarser  [4],  The  resulting 
temperature  distribution  obtained  from  the  corresponding  analytical  thermal  solution  is  continuous,  with 
the  only  approximation  being  the  averaging  of  the  power  densities  over  the  source-gate-drain  sublengths. 

The  fully  physical  LPM  makes  this  simple  prediction  of  device  performance,  Eq.  (5),  based  purely 
on  specification  of  structural  parameters,  such  as  semiconductor  layer  compositions,  widths  and  doping 
levels,  and  gate  recess  depth  and  width,  without  any  prior  experimental  characterisation.  The  large-signal 
RF  description  makes  a  similar  single  subroutine  call  to  the  LPM,  returning  instantaneous  drain-source 
voltage  as  a  function  of  instantaneous  values  of  source-gate  voltage  and  source  current.  The  main  physical 
difference  in  the  large-signal  RF  description  is  inclusion  of  a  displacement  contribution  to  the  current, 
which  provides  a  description  of  transistor  capacitances. 
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3.1.1  Low  computational  cost  of  the  LPM 


The  iterative  solution  of  the  electro-thermal  equation,  Eq.  (4),  requires  repeated  calls  to  the  LPM  via 
6i(Pi)  for  the  active  elements  i  =  1,  Despite  the  speed  of  the  quasi-2D  LPM  this  can  still  be  a 

significant  computational  expense.  However,  this  computational  cost  can  be  reduced  by  noting  that  many 
of  the  calls  to  the  LPM  can  occur  for  the  same  bias  point,  {VDS)  Vgs),  at  similar  values  of  the  channel 
temperature,  T. 

For  illustration,  assume  the  elements  i  to  be  identical,  e.g.  a  MMIC  grid  array  for  spatial  power 
combining  [54] ,  with  transistor  fingers  of  identical  layer  structure  and  design,  repeated  between  identical 
transistors  on  a  MMIC  and  between  identical  MMICs  in  the  array,  so  that  6i(Pi)  =  0(P<)  for  all  i. 
Then  for  a  given  bias,  the  LPM  can  be  called  for  a  range  of  power  dissipations,  Ph  sufficient  to  span 
the  range  of  Pi  generated  in  the  iterative  solution.  By  storing  the  resulting  values  6(Pi)  and  making  a 
simple  interpolation  between  them,  an  essentially  analytical  description  6(P)  can  be  generated  for  all 
P  of  interest,  at  the  cost  of  a  small  amount  of  memory  storage  for  each  bias  point.  This  approach 
avoids  wasteful  repeated  calls  to  the  LPM  in  the  vicinity  of  the  same  0* ,  Pi  during  iterative  solution. 
These  analytical  $(P)  can  be  calculated  prior  to  the  fully-coupled  electro-thermal  solution.  They  will  be 
generated  rapidly,  due  to  the  speed  of  the  LPM  and  provide  the  fully-coupled  electro- thermal  solution 
with  all  the  physical  predictive  power  and  accuracy  of  the  LPM,  at  the  minimal  computational  cost 
of  an  essentially  analytical  function  evaluation  and  small  memory  storage.  The  9(P)  will  also  provide 
essentially  analytical  derivatives  for  use  in  Newton-Raphson  solution  of  the  coupled  non  linear  equation, 
Eq.  (4). 

3.2  The  global  thermal  description 

As  the  electrical  part  of  this  coupled  model  is  fully  described  by  the  existing  LPM,  the  rest  of  the 
paper  concentrates  on  the  description,  construction  and  demonstration  of  the  global  thermal  model.  Key 
considerations  in  the  construction  of  the  global  thermal  model  were  as  follows: 

1.  It  must  provide  rapid  thermal  updates  in  coupled  electro-thermal  solutions  whilst  providing  a 
physical  (hence  predictive)  and  accurate  description  of  steady-state  heat  flow. 

2.  It  must  require  no  core  changes  to  the  LPM. 


10 


3.  It  must  treat  the  non  linearity  due  to  the  temperature  dependence  of  material  parameters. 

4.  It  must  allow  inclusion  of  radiative  and  convective  surface  fluxes  in  the  description  of  large  area 
systems. 

5.  It  must  allow  treatment  of  totally  arbitrary  layouts  and  layer  structures,  providing  a  flexible  tool 
for  electro-thermal  CAD. 

To  this  end  the  analytical  thermal  resistance  matrix  approach  was  adopted.  This  approach  is  described 
by  the  key  equation 

A  Oi  =  ^  RrHij  Pj  (7) 

3 

which  states  that  temperature  rise,  A0*,  at  active  device  element  i,  is  determined  entirely  by  the  power 
dissipations  Pj  at  all  active  device  elements,  j  =  1,  with  matrix  Q TR  determined  prior  to  the 

coupled  electro-thermal  calculation  purely  from  structural  information.  The  construction  and  demonstra¬ 
tion  of  the  consequences  of  such  a  relation,  Eq.  (7),  for  arbitrary  MMIC  and  hybrid  MIC  geometries,  form 
the  subject  of  the  rest  of  this  paper.  Key  features  of  the  analytical  thermal  resistance  matrix  approach 
are  indicated  below. 

Firstly,  as  thermal  updates  in  this  description  reduce  to  small  matrix  multiplications,  the  coupled 
electro-thermal  solution  is  rapid.  The  thermal  resistance  matrix  is  obtained  from  analytical  solutions  of 
the  heat  diffusion  equation  with  a  linear  ‘radiation5  boundary  condition,  so  it  treats  both  steady-state  and 
time-dependent  cases  with  surface  radiation  and  convection.  Non  linearity  due  to  temperature  dependent 
thermal  conductivity  is  treated  by  means  of  the  Kirchhoff  transformation  [55].  Arbitrary  3-dimensional 
geometries  are  treated  by  interface  matching  of  analytical  solutions  for  thermal  subelements  (such  as 
MMIC  dies,  air  bridge  legs  and  top,  flip  chips  and  solder  bumps,  and  other  surface  metallisation). 

Use  of  the  method  means  that  no  redundant  temperature  information  is  generated  within  the 
body  of  the  device.  Temperatures  are  only  obtained  at  the  active  devices,  as  required  for  the  coupled 
electro- thermal  solution.  (Though,  as  thermal  solutions  are  analytical,  once  active  device  power  dissipa¬ 
tions  have  been  obtained  self-consistently,  temperatures  can  be  obtained  at  any  point  within  the  volume 
to  any  degree  of  accuracy.  This  is  of  value,  for  instance,  in  model  validation  by  thermal  imaging  of  surface 
temperatures.) 

No  surface  or  volume  mesh  is  required  in  the  analytical  RTH  method.  Active  elements  and  thermal 
subsystem  interface  discretisations  are  specified  by  the  coordinates  of  the  elementary  area  boundaries. 
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The  thermal  resistance  matrix  description  R  th  gives  a  one-to-one  equivalent  circuit  model  for  the 
thermal  interactions  [17],  Figure  4,  which  enables  the  immediate  incorporation  of  the  coupled  electro¬ 
thermal  description  into  general  electrical  network  based  solvers. 

The  form  of  the  equivalent  circuit  illustrated  above  [17,  18]  lends  itself  readily  to  extension  to  arbi¬ 
trarily  large  numbers  of  nodes.  This  contrasts  with  some  more  conventional  thermal  equivalent  circuit 
representations  which  attempt  to  identify  directly  circuit  elements  with  obvious  physical  parameters,  and 
which  can  grow  rapidly  more  convoluted  with  increase  in  size  [25].  The  equivalent  circuit  description  illus¬ 
trated  above,  and  most  immediately  identified  with  the  thermal  resistance  matrix  formulation  presented 
here,  lends  itself  readily  to  circuit  description  of  complex  thermal  problems. 

4  Construction  of  the  Resistance  Matrix 

Having  presented  the  motivation  for  the  analytical  thermal  resistance  matrix  approach,  Eq.  (7)  is  now 
derived  for  the  thermal  steady-state  case  and  RTH  obtained  for  the  basic  example  of  a  homogeneous, 
heatsink  mounted  MMIC,  bearing  an  arbitrary  distribution  of  power  transistor  fingers. 

This  is  followed  by  discussion  of  RTH  for  a  MMIC  comprised  of  multi-layers  of  differing  thermal 
conductivity,  which  includes  the  cases  of  uniform  surface  metallisation,  metal  coated  MMIC  base,  and 
doping  dependent  thermal  conductivities  in  MMIC  layers. 

Treatment  of  MMIC  superstructure  is  then  presented  followed  by  discussion  of  the  construction  of 
R  using  a  double  Fourier  series  finite  element  method  for  irregular  surface  metallisation.  This  is 
followed  by  construction  of  RTH  for  an  inhomogeneous  MMIC  with  vias  or  partial  substrate  thinning. 

4.1  Homogeneous  MMIC 

The  steady-state  heat  diffusion  equation,  in  the  absence  of  volume  heat  sources  or  sinks,  is  given  by 

V.  [k(T)VT]  =  0,  (8) 

where  T  is  physical  temperature  and  k{T)  is  temperature  dependent  thermal  conductivity.  This  equation 
is  non  linear  through  n{T).  The  equation  can  be  linearised  by  the  Kirchhoff  transformation  [55] 

6  =  Ts  +  —  fT  n(T)dT,  (9) 

KS  JTs 
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where  Ks  =  k(Ts).  This  gives  the  linearised  steady-state  heat  diffusion  equation 


V26  =  0. 


Then  obtaining  the  analytical  solution  for  linearised  temperature  9 ,  physical  temperature  T  is  imme¬ 
diately  given  by  T  =  T(9 ),  where  function  T(0)  is  determined  from  the  Kirchhoff  transformation,  Eq. 


(9).  For  k(T)  of  the  form 


/  rp  \  — b 

k(T)  =  KS  (  jT  )  > 


appropriate  to  GaAs  in  temperature  ranges  of  interest  [56],  T(9)  is  given  by  [26] 


T(0)  =  Ts  1- 


(b-l)(9-Ts) 


To  generate  a  solution  of  Eq.(10),  the  physical  system  must  be  specified  and  surface  boundary  condi¬ 
tions  must  be  imposed.  Eq.  (10)  can  be  solved  for  general  linear  boundary  conditions,  but  for  simplicity 
adiabatic  boundary  conditions  are  assumed  on  the  side  faces  of  the  MMIC.  Active  device  power  dissipa¬ 
tion  and  heatsink  mounting  are  described  by  a  generalised  ‘radiation’  boundary  condition  on  the  MMIC 
top  and  bottom  surfaces.  For  illustration  of  the  coordinate  system,  the  structure  of  Figure  5  is  considered. 
This  consists  of  a  GaAs  MMIC,  0  <  x  <  L,0  <  y  <  W90  <  z  <  D,  bearing  an  arbitrary  number  of 
surface  heating  elements,  areas  A,  dissipating  (surface  average)  power  densities,  Pi  (corresponding  to 
both  active  and  passive  devices). 

Solving  the  linearised  steady-state  heat  diffusion  equation  subject  to  adiabatic  flux  boundary  condi¬ 
tions  on  the  MMIC  side  faces,  and  the  generalised  linear  radiation  boundary  condition, 

d0 

+  H0,d  [0  -  0O)D(x, y)}  +po ,d(x’V)  =  °>  (13) 

on  the  MMIC  top  and  bottom  surfaces,  the  full  solution  for  the  linearised  temperature  by  separation  of 


variables  is 


6  =  A  +  Bz  4-  cos  Xmx  cos  \iny 


X  {pmn  COSh  'yrnnZ  +  Smn  sinh7mnz)> 


where  Yl!mn  means  sum  over  m,  n  =  0, 1, 2, ...  excluding  (m,  n)  =  (0,  0), 


13 


(16) 


+  /4> 


1  rL  rw 

=  ~Jjw  J  J  Pd{x,v)  -  HD9D(x,y)dxdy 


-((Xdks  +  HdD)B, 


Ksa0B 


1  fL  fw 

WV  J  J  MX,V)  ~  Ho8o(x,y)dxdy 


- H0A , 


-'mn  —  'Tmn  ^mn^Q  ^5 

Si  fpV  COS  \mx  cos  p„y  [p0 (x,  y)  -  ff0flo (a:,  y)]  dxdy 
^-(1  +  <5mo)(l  +  fino) 


Cmn  [aDKS'rmn  sinh7mnD  +  HD  cosh  7 mnD) 

"h  Smn  [aDKslmn  cosh  7 m„C  +  HD  sinh7mnZ>] 

fL  fw  cos  Xm.x  cos ^nylPD(x,y)—Hr)OD(oc,y)]dxdy 
=  ~  _a_a  ^(l+^oKl+ino)  ‘ 

Choosing  ao  =  1,  otp  =  0,  H0  =  0,  Ho  =  1,  Po(x,y)  =  Y^i^i(x>y)pi>  Pd{x,v)  =  0,  6o{x,y)  =  0, 
0D(x  ,y)  =Ts  describes  a  homogeneous,  heatsink  mounted  MMIC  at  base  temperature  Ts  with  imposed 
surface  flux  Po(x,y)  due  to  active  device  elements,  and  no  surface  convective  or  radiative  losses.  Here 
Si(x,  y)  equals  1  in  elementary  areas  Di  dissipating  uniform  power  densities  P*,  and  0  otherwise. 

This  gives  the  explicit  form  for  the  linearised  temperature  throughout  the  MMIC, 

9  -  ’■-SfEw* 

/ 

5>  Am#  cos  pny 

mn 

x(sinh7mnz  -  tanh7mnjDcosh7mn^) 


KsLW  (1  +  £mo)(l  +  &no)lv 


Ep  p 


where  elementary  area  integral  Pmn  is  defined  by 


i„  =  J  J  cos  Xmx cos yny  dxdy 
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and  the  standard  result, 


L 

J  COS  Xmx  COS  A mfX  dx  =  ~Smmf  (1  +  5 


(24) 


has  been  used  (and  similarly  for  the  corresponding  y  integral)  with  6mn  the  Kronecker  delta  function. 
To  construct  the  thermal  resistance  matrix,  surface  averaged  temperature  rises,  A Biy  are  evaluated  as 


A  0i  = 

Then  Eq.  (7)  is  obtained  immediately  with 

RrHii  =  ksLW  [' 


_  /  JDi  Q  —  Ts  dxdy 


I  JDidxdy 


(25) 


Dll  o  + 


E 


4  tanh  7mn  D  PmnI3m, 


7mn  (1  +  ^mo)(l  “h  <^no)  7q 


00 


(26) 


This  expression  illustrates  the  main  points  of  the  analytical  thermal  resistance  matrix  approach.  The 
thermal  resistance  matrix  is  evaluated  just  once,  prior  to  coupled  electro-thermal  calculation,  purely  in 
terms  of  structural  parameters.  As  it  is  obtained  for  linearised  temperature,  it  is  independent  of  both 
temperature  and  power,  and  hence  independent  of  electrical  bias  point.  With  the  non  linear  relation 

Mi  =  AQi(Pi),  (27) 

from  the  physical  active  device  model,  combination  with  the  global  thermal  description,  Eq.  (7),  gives 

A  ei(Pi)  =  '£RTHijPj  (28) 

3 

which  is  a  small,  simple,  non  linear  system  to  be  solved  self-consistently  for  the  power  densities,  Pj. 
Having  obtained  the  Pj  from  solution  of  the  coupled  electro-thermal  problem,  for  instance  by  Newton 
methods,  the  full  electrical  solution  is  obtained.  Also,  the  temperature  at  any  point  within  the  MMIC  is 
then  given,  if  required,  by  Eq.  (22). 

As  an  illustration,  the  surface  temperature  distribution  calculated  using  the  coupled  electro-thermal 
HEMT  model  is  shown  in  Figure  6  for  a  power  HEMT  with  60  fingers.  Total  power  dissipated  is  2.8  W, 
Ts  =  300  K,  and  plot  temperature  varies  from  below  36  °C  to  53  °C.  Calculation  of  a  full  set  of  I-V 
curves  for  this  device  took  some  10’s  of  minutes  based  on  a  non  optimised  relaxation  algorithm. 

The  simple  expression,  Eq  (26),  represents  a  full  3-dimensional,  solution  of  the  heat  diffusion  equation 
for  a  MMIC  with  an  arbitrary  surface  distribution  of  power  transistors  (and  heat  dissipating  passive 
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elements),  each  power  transistor  having  an  arbitrary  grouping  and  arrangement  of  transistor  fingers.  It 
includes  fully,  finite  length  and  end  effects  for  heat  diffusion  from  the  active  layers,  and  treats  the  finite 
volume  effects  of  the  die  without  any  need  for  the  simplifying  assumption  of  an  infinite  or  semi-infinite 
substrate.  The  extension  to  other  realisations  of  the  radiation  boundary  condition,  Eq.  (13),  for  instance 
deriving  simple  resistance  matrix  expressions  for  large  area  substrates  with  convective  and  radiative 
surface  losses,  is  immediate. 

Use  of  the  expression,  Eq.  (26),  in  coupled  electro-thermal  simulations  of  power  FETs  and  MMICs, 
neglects  the  effects  of  surface  metallisation  and  resulting  heat  spreading.  It  will  therefore  generally 
overestimate  the  magnitude  of  KTH >  producing  simulated  device  temperatures  higher  than  those  realised 
experimentally.  Simple  use  of  Eq.  (26)  can,  therefore,  only  to  be  expected  to  produce  a  safe  upper  limit  on 
temperature  rises.  It  will  give  correct  order  of  magnitude  estimates  for  temperature  rises,  and  demonstrate 
correct  trends  in  temperature  variation  with  changes  in  parameters,  but  accurate  prediction  of  RTH 
based  purely  on  the  physical  thermal  model  will  typically  require  description  of  metallic  superstructure, 
as  described  in  the  sections  following. 


4.2  Generalisation  of  RTH 


In  Eq.  (25)  of  the  previous  section,  surface  average  temperature  rises  were  calculated  for  the  elementary 
areas,  D{,  over  which  the  active  device  power  densities,  P*,  were  dissipated.  As  temperature  rises  in  a 
multi-finger  device  will  typically  be  highest  in  the  immediate  vicinities  of  peak  power  dissipations,  the 
thermal  resistance  matrix  of  Eq.  (26)  will  return  peak  temperatures. 

For  use  in  the  coupled  electro-thermal  formulation,  with  the  LPM  employing  a  uniform  channel 
temperature  approximation,  the  active  device  temperature  that  is  fed  into  the  electrical  model  should 
be  the  average  temperature  over  the  whole  source-gate-drain  area.  This  means  that  the  construction  of 
surface  average  temperature  rises,  as  in  Eq.  (25),  should  be  modified  so  that  temperatures  are  averaged 
over  source-gate-drain  areas,  D corresponding  to  each  active  device  power  dissipating  area,  Di. 

This  gives  immediately  for  the  generalised  thermal  resistance  matrix, 


RtHu 


^K+ 

4  tanh7mnP 

~f  7mn(l  +  ^mo)(l  +  <5no) 


Ifi  p 

-Lmnmn 


(29) 
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where  elementary  area  integrals  J*n,  /qo  are  still  of  the  form^  Ech  (23)>  but  now  integrations  are  over 
areas  D[. 

More  generally,  having  obtained  an  analytical  solution  of  the  form,  Eq.  (22),  temperature  rises  can 
be  calculated  at  any  point,  or  averaged  over  any  area  or  volume,  irrespective  of  any  correlation  to  power 
dissipating  areas.  Thermal  resistance  matrices  so  constructed  will  generally  be  of  a  different  explicit  form 
to  Eqns.  (26)  and  (29)  and  will  not  generally  be  square. 

Use  of  source-gate-drain  areas,  D in  the  thermal  resistance  matrix,  is  to  be  distinguished  from  use  of 
an  effective  ‘hot  strip’  width  much  wider  than  the  dimensions  of  the  physical  power  dissipating  area,  as  in 
[16].  In  [16],  a  wide  active  power  dissipating  area  was  assumed  to  obtain  agreement  with  low  resolution 
thermal  measurements.  In  the  formulation  presented  here,  the  extents  of  the  power  dissipating  regions, 
Di ,  are  physical  and  are  obtained  from  the  LPM.  The  use  of  source-gate- drain  averaged  temperature 
rises  over  areas,  D[ ,  is  required  for  physical  consistency  with  the  uniform  channel  temperature  version 
of  the  LPM.  The  effects  of  low  resolution  in  thermal  measurements  can  be  simulated  by  averaging  the 
analytical  result,  Eq.  (22),  over  pixel  sized  areas,  after  generation  of  the  self-consistent  electro- thermal 
solution. 

4,3  Multilayer  MMIC 

The  simple  description  of  the  homogeneous  MMIC,  presented  above,  is  readily  generalised  to  treat  multi¬ 
layer  systems  by  use  of  a  transfer  matrix,  or  two-port  network,  approach  [9].  This  is  based  on  matching 
of  Fourier  components  at  interfaces,  and  corresponds  to  use  of  the  double  cosine  transform  to  convert  the 
3-dimensional  partial  differential  equation,  Eq.  (10),  into  a  1-dimensional  ordinary  differential  equation 
for  the  ^-dependent  double  Fourier  series  coefficients.  Matching  of  linearised  temperature  and  flux  at 
the  interfaces  of  a  multi-layer  structure  can  then  be  imposed  by  use  of  a  2  x  2  transfer  matrix  on  the 
Fourier  series  coefficients  and  their  derivatives.  Arbitrary  N-level  structures  can  be  treated.  Different 
thermal  conductivities  can  be  assumed  in  each  layer  allowing  treatment  of  composites  like  Cu  on  AIN 
(both  having  temperature  independent  thermal  conductivities)  and  MMICs  with  conductivities  varying 
from  layer  to  layer  due  to  differences  in  doping  levels  (all  layers  having  the  same  functional  form  for 
the  temperature  dependence  of  the  conductivity).  The  method  can  be  generalised  further,  by  imposing 
specified  flux  discontinuities  at  the  interfaces.  The  solution  then  represents,  for  instance,  the  case  of  a 
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MMIC  with  active  device  channel  buried  by  a  thin  layer  of  semiconductor. 

For  instance,  consider  a  2-layer  structure  of  the  general  form  of  Figure  5,  but  with  top  layer  0  < 
z  <  Du  thermal  conductivity  «i  and  bottom  layer  <  z  <  D2,  thermal  conductivity  «2.  With 
power  dissipations  Pi  over  elementary  areas  Di  at  the  interface,  z  =  D±,  the  (linearised)  temperature 

distributions  6^  and  in  the  top  and  bottom  layers  respectively,  are  given  by 

/ 

on)  =  Ai  d+£  COS  \mX  COS  fJLny  ^mn  cosl1 7 rnn^j  (30) 

mn 

0(2)  =  Ts  +  B(2)  [z  -  (D1  +D2)}  + 

/ 

E«  \m,X  cos  finyS%l  X 

ran 

[sinh7mnz  -  tanh 7mn (D\  +  D2)  cosh7mnz] , 


where, 


A™  =  Ts-B^Lh, 

^mn  ~  mn  [tanll7mn.Di  ~  tanh 7mn (D\  +  D2)\  ■ 


B(2)  = 


s£>  = 


4  cosech7mn-Di 
LW  (1  +  SmQ  )  ( 1  +  Sno)jn 


£ ^n-Pi 


iti  [tanh7m„Di  -  tanh7m„(Z>i  +  £>2)] 
-«2  [cotli7mnDi  -  tanh7mn(Di  +  D2)] 


The  expression  for  the  corresponding  thermal  resistance  matrix,  describing  the  temperature  rises  of  the 
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power  dissipating  elements  at  the  interface,  is 


RTHij  =  ^ZW  [D2I00+ 

_ 4 _  IlnnH, 

2-^mn  (l+(5Tno)(H“<5no)'ymn  Iqq 


(36) 


tanh7mrt  (D1+D2)— tanh  7mn£>i 


1  -  tanh2  jmnDi 


“  (l  -  «)  tanh7mnpi  +  D2)  tanh7mnDi  j 


Adopting  an  effective  value  for  «i,  to  allow  for  the  different  functional  forms  of  temperature  dependent 
thermal  conductivity  in  metal  and  GaAs,  this  expression  can  be  used  to  provide  a  simple  approximation 
to  the  thermal  resistance  matrix  of  a  heavily  metallised  power  FET  or  MMIC. 


4.4  MMIC  superstructure 

It  has  been  demonstrated  that  inclusion  of  surface  metallisation  is  essential  for  accurate  description  of 
thermal  effects  in  power  devices  [4,  5,  7],  Comparison  with  experiment  for  multi-finger  power  HBTs 
shows  that  the  simple  thermal  description  corresponding  to  the  resistance  matrix  of  Eq.  (26)  is  highly 
accurate  when  combined  with  a  simple  model  of  heat  shunting  by  an  air  bridge  [17].  The  extension  of 
the  analytical  thermal  resistance  matrix  approach  to  include  descriptions  of  surface  metallisation  and  air 
bridges,  and  other  vertical  geometries  such  as  flip  chips  and  solder  bumps,  is  now  presented. 

The  extension  to  include  surface  superstructure  is  achieved  by  solving  the  heat  diffusion  equation 
analytically  for  thermal  sub  elements,  then  combining  thermal  resistance  matrices  for  subsystems  by 
matching  of  temperature  and  flux  at  discretised  interfaces.  In  the  case  discussed  below,  the  multi-layer 
problem  is  solved  for  discretised  temperatures  at  thermal  sub  element  interfaces.  However,  the  matching 
problem  can  equally  be  cast  in  terms  of  solution  for  discretised  interface  fluxes.  The  former  method 
is  described  here,  as  the  latter  produces  a  more  complicated  formulation  for  the  system  of  equations 
to  be  solved  for  the  interface  unknowns.  (This  is  because  specification  of  flux  boundary  conditions  on 
top  and  bottom  surfaces  of  a  thermal  sub  element  only  defines  temperature  upto  an  arbitrary  constant, 
and  imposes  a  subsidiary  energy  conservation  constraint  on  total  flux  into  and  out  of  the  volume.  When 
matching  sub  elements  at  interfaces,  these  energy  conservation  constraints  serve  to  determine  the  arbitrary 
constants  in  the  temperatures,  but  require  moderately  complicated  manipulation  of  the  simultaneous 
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equations  for  the  interface  unknowns,  to  cast  them  in  terms  of  only  independent  fluxes  and  the  related 
arbitrary  temperature  constants.) 

To  generate  analytical  solutions  allowing  interface  matching  of  thermal  sub  elements,  ofo,  otD)  #o> 
HD:  Po(x,y),  pD{x,y)i  80(xyy)  and  6o(x,y)  of  Eq.  (13)  are  specified  as  in  Section  4.1,  but  withp0(^,2/) 
written  Po(x,y)  =  YliSoi(x,y)Poi  and  the  uniform  MMIC  base  temperature  replaced  by  the  discretised 
temperature  distribution  #£>(£,  y )  =  JE  &Dj  ix>  y)@Dj*  Similarly  for  thermal  sub  element  surface  temper¬ 
ature,  8oi,  and  base  flux,  Poj-  Then  following  the  procedure  of  Section  4.1,  specifying  either  discretised 
base  temperature  or  base  flux,  the  following  relations  are  obtained: 


Ooi  =  X/  R THu '  p0i'  +  X/  ^$£>3  > 

i '  j 

(37) 

PDj  =  E  ^THjiPoi  +  ^2  Tjj’6Dj’, 

(38) 

where  Rth ti,  is  given  by  Eq.  (26)  and 
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(39) 

/°jn  and  are  the  area  integrals  of  Eq.  (23)  evaluated  over  elementary  areas  D*  on  the  surface,  z  =  0, 
and  over  elementary  areas  Dj  on  the  base,  z  ~  D,  respectively.  N  is  the  number  of  discretised  base 
elementary  areas,  Dj ,  and  a  takes  any  value  with  dimensions  W/(K.m2). 

To  illustrate  the  interface  matching  approach,  the  global  thermal  resistance  matrix  is  now  constructed 
for  the  case  of  N  pieces  of  rectangular,  but  otherwise  arbitrary,  metallisation  on  the  surface  of  an  otherwise 


20 


homogeneous  heatsink  mounted  MMIC.  For  simplicity  of  illustration,  only  a  single  metal  layer  is  consid¬ 
ered  (though  each  individual  piece  of  surface  metal  can  be  of  any  required  thickness).  The  generalisation 
of  the  this  approach  to  multilayer  structures  such  as  air  bridges  and  flip  chips  is  immediate. 

For  a  single  metal  layer,  surface  flux  terms,  Poi>  in  Eqns.  (37)  and  (38)  are  identically  zero.  As  the 
surface  temperature  of  the  metal,  60i ,  is  not  required  for  the  coupled  electro-thermal  solution,  only  Eq. 
(38)  is  required.  Written  in  matrix  notation  this  becomes  simply 


pfa)  _  p(n)^(n) 


(40) 


for  each  piece  of  surface  metallisation,  n,  with  n  =  1  ,  ...,iV.  The  thermal  conductivity  ks  in  Z(n)  >  Eq.  (39), 
for  each  piece  of  metallisation,  n,  can  be  chosen  at  a  different  representative  Kirchhoff  transformation 
temperature,  Ts-  Such  a  choice  allows  more  accurate  treatment  of  the  difference  in  functional  form 
between  conductivity  in  GaAs  and  metal  [57].  This  technique  is  described  in  more  detail  in  Section  4.5 
for  the  inhomogeneous  MMIC  with  vias. 


In  matrix  notation,  Eq.  (7)  for  the  MMIC  die  becomes 
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es-Tss 

6*  -  Tsi 


nss  jysi 

£th  =th 
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V  =TH  =TH  /  V  -  / 
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(  o,  \ 
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(41) 


where  Tss  and  Tsf  are  constant  vectors  with  all  elements  equal  to  Ts,  and  the  matrix  equation  has  been 
partitioned  by  active  device  elementary  surface  areas,  s,  and  interface  elementary  areas  between  MMIC 
die  and  metal,  i. 

Matching  flux  and  (linearised)  temperature  at  the  interface  between  metal  and  MMIC  die,  the  fol¬ 
lowing  relation  is  obtained 

Ms  =RS^IP\  (42) 


where  Ps  is  the  vector  of  MMIC  active  device  power  dissipations,  R910^  is  the  global  thermal  resistance 
matrix  for  the  coupled  GaAs  and  metal  system,  and  A6S  is  the  vector  of  MMIC  active  device  temperature 
rises. 

The  global  resistance  matrix  is  given  explicitly  by 


R9i°Z  =  Rs°  +R%„RRi‘ 


(43) 


where 

E  =  diag(£(1\...,l(n),-..,lW) 
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X  [r  -  i”Hdiag  ...,T(jV))] 


(44) 


P  is  the  identity  matrix,  and  temperature  rises  A0S  are  now  defined  by 

A 6s  =  0s  -  (ls  +  EStHE)  Is*  ,  (45) 

where,  again,  Is  is  the  identity  matrix. 

Eq.  (42)  illustrates  the  main  points  of  the  hierarchical  construction  of  global  thermal  resistance 
matrices  for  complex  systems.  Firstly,  the  order  of  the  global  thermal  resistance  matrix,  R9^,  is  small, 
determined  only  by  the  number  of  active  device  elements,  not  the  level  of  internal  discretisation.  This 
implies  rapid  temperature  updates  in  the  iterative  coupled  electro-thermal  solution.  Secondly,  R9^  is 
constructed  purely  from  simple  matrix  manipulations  on  elementary  matrices  given  by  simple  analytical 
expressions  for  thermal  subsystems.  The  global  resistance  matrix  is  constructed  just  once,  prior  to 
the  coupled  electro-thermal  solution,  so  its  construction  has  no  impact  on  coupled  run  time.  The  only 
significant  approximation  in  construction  of  R9^  is  the  level  of  interface  discretisation.  This  is  determined 
only  by  available  computing  power  and  required  accuracy,  and  description  by  R9^  becomes  exact  in  the 
limit  of  infinitely  fine  interface  discretisation.  Bonani  et  al  have  typically  found  1000  -  2000  nodes  to  be 
sufficient  in  their  hybrid  finite  element  Green’s  function  description  of  MMIC  surface  metallisation  [2,  3]. 

4.4.1  Non  linear  interface  matching 

In  the  hierarchical  treatment  of  MMIC  superstructure  just  presented,  linearised  temperatures  were 
matched  at  subsystem  interfaces.  This  is  only  strictly  justified  when  the  functional  form  of  the  temper¬ 
ature  dependence  of  the  thermal  conductivity  is  the  same  in  both  matched  layers.  When  the  functional 
forms  are  different,  matching  of  physical  temperature,  T,  imposes  a  non  linear  relation  on  the  linearised 
temperatures,  6 ,  in  the  two  subsystems.  This  arises  as  a  result  of  the  different  relations  T  =  T(0)  obtained 
from  the  respective  Kirchhoff  transformations,  Eqns.  (9),  (11)  and  (12).  For  GaAs,  exponent  b  in  Eq. 
(11)  is  typically  equal  to  1.22  [56]  whereas  for  metal  b  equals  0,  i.e.  the  thermal  conductivity  of  metal  is 
essentially  temperature  independent.  This  difference  means  that  matching  of  linearised  temperatures  at 
the  GaAs/metal  interfaces  implies  a  non  linear  matching  relation.  (Matching  of  fluxes,  however,  remains 
linear  after  differing  Kirchhoff  transformations.) 
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This  non  linear  interface  matching  can  be  avoided  to  a  good  approximation,  by  assuming  the  same 
functional  form  for  the  temperature  dependent  thermal  conductivity  of  the  metal  and  of  the  GaAs.  The 
resulting  error  in  surface  temperatures  is  calculated  to  be  typically  a  few  per  cent  [7,  57].  This  error  can  be 
reduced  further  in  multi-component  metal  on  GaAs  systems,  by  assuming  a  different  Kirchhoff  transfor¬ 
mation  temperature,  Ts,  for  each  piece  of  surface  metallisation.  Choosing  these  Kirchhoff  transformation 
temperatures  to  give  a  thermal  conductivity  &s  close  to  the  metal  thermal  conductivity,  minimises  the 
errors  introduced  by  adoption  of  an  incorrect  functional  form  for  the  temperature  dependence. 

In  cases  where  the  non  linear  interface  matching  cannot  be  neglected,  the  thermal  resistance  matrix 
approach  allows  formulation  of  a  non  linear  system  of  equations  for  the  correctly  matched  temperatures 
[19].  For  the  case  of  a  grid  array  of  identical  GaAs  MMICs  on  an  AIN  substrate  (with  essentially 
temperature  independent  thermal  conductivity),  arguments  similar  to  those  of  Section  4.4  give  rise  to 
the  non  linear  system  for  interface  temperatures  0^  of  MMICs  n  =  1, ...,  M, 

NL(6 &>)  =  J2&THnm  (— TJf  ^  +  U(D  ] ')  ,  (46) 

m 

in  terms  of  the  MMIC  surface  flux  densities  P  NL(0)  represents  the  non  linear  relation  between  tem¬ 
peratures  matched  at  the  interfaces.  The  substrate  thermal  resistance  matrix  has  been  partitioned 
into  submatrices  R ,  with  n,  m  =  1, M  running  over  the  individual  MMICs,  and  S  and  T  are 
defined  as  in  Eqns.  (39).  Solving  self-consistently  for  the  interface  temperatures,  Eq.  (37)  provides  a  non 
linear  relation  between  surface  temperatures  6  ^  and  surfaces  fluxes  P  ^  for  the  MMIC  grid  array. 

This  non  linear  problem  for  the  interface  temperatures  can  again  be  solved  by  Newton  methods. 
The  number  of  unknowns  is  determined  by  the  level  of  discretisation  of  the  MMIC  bases,  and  can  be 
minimised  by  choosing  the  shapes  of  the  elementary  areas  to  reflect  the  expected  behaviour  of  the  MMIC 
base  temperature  profiles. 

Figures  7  and  8  illustrate  the  simulation  of  a  5  x  5  power  FET  array,  with  explicit  non  linear  matching 
of  MMIC/substrate  interface  temperatures.  The  MMIC  array  substrate  was  heatsink  mounted  and 
adiabatic  surface  boundary  conditions  were  assumed.  Single  surface  averaged  MMIC  base  temperatures 
were  matched  at  the  MMIC-substrate  interfaces.  The  coupled  electro-thermal  simulation  assumed  non 
linear  matching  of  linearised  temperatures  at  the  interfaces,  due  to  the  different  functional  forms  of  the 
temperature  dependences  in  GaAs  MMICs  and  AIN  substrate.  Coupled  electro-thermal  simulation  time 
for  1  DC  bias  point  was  ~10  minutes  (after  construction  of  the  respective  thermal  resistance  matrices). 
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This  was  based  on  a  simple  non  optimised  relaxation  algorithm  and  repeated  calls  to  the  LPM  in  the 
iterative  coupled  solution.  Use  of  Newton  methods,  combined  with  prior  calculation  and  interpolation  for 
the  LPM  temperature  dependence,  is  expected  to  produce  significant  further  speed-up  up  in  simulation 
run  time  which  was  dominated  by  calls  to  the  LPM. 

4.4,2  Non  linear  surface  fluxes 

Surface  radiative  fluxes  are  unimportant  for  the  small  areas  of  heatsink  mounted  MMICs  and  even  hybrid 
MICs.  Convective  fluxes  are  non  trivial  to  calculate  for  small  areas  and  fine  surface  structure.  However, 
for  large  area  substrates,  from  MMIC  grid  arrays  upto  circuit  board  level,  convective  and  radiative  losses 
are  known  to  become  significant  and  at  high  powers  and  temperatures  surface  flux  non  linearities  need 
to  be  considered.  The  linear  thermal  resistance  matrix  approach  allows  treatment  of  non  linear  surface 
fluxes  by  discretising  the  whole  substrate  surface  (not  just  the  active  device  elements)  and  treating  the 
radiative  and  convective  surface  fluxes  as  unknown  imposed  fluxes  to  be  obtained  self-consistently. 

Adopting  this  approach,  Eq.  (7)  gives  immediately  the  following  non  linear  relation  for  the  surface 
temperature  rises  A $i  [19], 

A0*  =  ^2  Rth^  Pj  ~  Y2  QTHir  nK®r)  - 

3  r 

(47) 

s 

in  terms  of  the  imposed  surface  fluxes  Pj.  nl(6)  is  the  surface  flux  non  linearity,  0r  the  substrate  surface 
temperature  on  the  front  surface,  and  6rsev  substrate  temperature  on  the  reverse  surface.  Qth ir  and 
VrHis  are  thermal  resistance  matrices. 

This  relation  leads  immediately  to  an  iterative  solution  for  treatment  of  surface  flux  non  linearity.  As 
the  non  linear  contribution  to  the  surface  flux  is  obtained  as  a  perturbation  about  the  linear  radiation 
boundary  condition,  and  as  the  linear  case  can  be  solved  to  provide  a  good  first  guess  in  the  iterative 
non  linear  solution,  the  problem  is  expected  to  be  rapidly  convergent  and  stable. 

The  required  level  of  surface  discretisation  and  hence  the  number  of  unknowns,  A 0iy  in  the  non  linear 
system,  will  depend  on  the  magnitude  of  the  substrate  thermal  conductivity  and  the  corresponding  level 
of  uniformity  of  the  substrate  temperature. 

Where  the  non  linearity  is  too  strong  to  allow  rapid  solution  via  the  iteration  of  the  linear  thermal 
resistance  matrix,  direct  solution  of  the  heat  diffusion  equation  for  arbitrary  non  linear  boundary  condi- 
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tions  can  be  achieved  by  a  combination  of  the  double  Fourier  series  expansion  solution,  Eq.  (14),  and  a 
collocation  technique.  This  solution  is  exact  to  the  extent  that  the  double  Fourier  series  has  converged. 
This  approach  produces  a  non  linear  thermal  resistance  function  Rth(P u  R, Fn)-  The  solution  is 
computationally  much  more  expensive  than  simple  linear  matrix  construction,  but  possibly  of  value  for 
high  accuracy  at  high  power  dissipations  and  temperatures.  The  linear  thermal  resistance  matrix  solution 
requires  surface  discretisation  and  corresponding  thermal  terminals  in  a  thermal  network  solution  for  non 
linear  surface  fluxes.  By  contrast,  the  non  linear  thermal  resistance  function  provides  a  description  of 
arbitrary  surface  flux  non  linearity  without  introducing  additional  surface  terminals  in  the  non  linear 
thermal  network. 

4.4.3  #___  from  finite  element  or  finite  difference  schemes 

—  I  H 

The  generation  of  the  thermal  resistance  matrix  has  so  far  been  described  on  the  basis  of  analytical 
solutions  of  the  heat  diffusion  equation.  However,  RTH  can  also  be  generated  experimentally  or  by 
numerical  simulation.  Direct  construction  of  RTH  based  on  physical  temperatures  obtained  from  non 
linear  experiment  or  simulation,  only  allows  linearisation  in  a  small  signal  sense,  and  is  of  very  limited 
applicability.  However,  if  the  system  non  linearity  can  be  attributed  almost  totally  to  temperature 
dependent  conductivity,  i.e.  if  surface  flux  non  linearities  are  negligible,  then  physical  temperatures 
can  be  converted  to  linearised  temperatures  via  the  Kirchhoff  transformation,  Eq.  (9),  and  a  thermal 
resistance  matrix  for  linearised  temperature  can  be  constructed.  The  thermal  description  so  obtained  is 
independent  of  temperature  and  power  dissipation,  and  so  of  device  bias  point.  It  is  not  restricted  to 
small  temperature  variations. 

Numerical  simulations  offer  an  increased  degree  of  flexibility  beyond  that  obtainable  by  experiment. 
Where  a  system  is  composed  of  subsystems  with  different  functional  forms  for  the  thermal  conductivity, 
each  subsystem  can  be  simulated  separately,  and  a  thermal  resistance  matrix  constructed  for  linearised 
temperature  in  each  case.  Solving  a  system  of  non  linear  equations  for  matched  interface  temperatures 
again  gives  a  thermal  description  valid  for  all  temperatures  and  power  densities  and  not  limited  to  small 
temperature  swings. 

To  generate  RTH  experimentally  or  numerically,  Eq.  (7)  gives  the  definition 


Generation  of  a  thermal  resistance  matrix  of  order  N  x  AT,  corresponding  to  N  active  device  elements, 
therefore  requires  N  individual  measurements  or  simulations,  dissipating  power  in  just  once  active  device 
at  a  time,  and  measuring  the  resulting  temperature  rises  at  all  active  devices.  (The  number  of  measure¬ 
ments  can  be  reduced  if  the  system  is  symmetrical.)  However,  to  generate  experimentally  or  numerically, 
a  thermal  resistance  matrix  giving  comparable  accuracy  to  a  matrix  obtained  from  analytical  expres¬ 
sions,  can  be  prohibitively  time  consuming.  For  example,  to  construct  RTH  for  a  power  MESFET  with 
10  fingers,  each  divided  into  active  device  elements  5  times  along  both  length  and  width,  would  require 
250  individual  simulations.  At  ~l/2  hour  per  simulation  for  a  finite  difference  calculation  [8],  the  time 
demands  are  heavy,  even  making  full  use  of  symmetry. 

4.4.4  Double  Fourier  series  finite  element  method 

To  treat  surface  metallisation  on  a  MMIC  die,  vertical  matching  of  rectangular  surface  elements  was 
described  in  Section  4.4.  However,  there  is  an  approximation  inherent  in  this  approach,  due  to  the 
imposition  of  adiabatic  side  wall  boundary  conditions  on  the  surface  elements.  These  conditions  mean 
that  no  direct  heat  flow  will  occur  between  adjoining  surface  metal  areas.  Generally  this  is  not  expected 
to  have  much  impact  on  calculated  thermal  resistance  matrices,  as  the  heat  flow  through  the  much  greater 
cross  sectional  area  of  the  MMIC  die  itself  is  unaffected  by  the  adiabatic  side  wall  boundary  conditions 
on  the  metal  elements.  Also,  any  adjoining  areas  of  metallisation  can  always  be  chosen  to  lie  in  low 
temperature  regions  away  from  centres  of  the  active  devices,  so  that  the  impact  on  heat  flow  of  artificial 
adiabatic  barriers  will  be  small. 

However,  at  the  expense  of  increased  algebraic  complexity,  the  analytical  thermal  resistance  matrix 
approach  immediately  generalises  to  remove  this  approximation.  The  thermal  description  of  a  cuboid 
with  arbitrary  fluxes  defined  on  all  6  faces  is  given  by  the  general  solution 

9  =  c0  +  cxx  +  cyy  +  czz  +  cxxx2  +  cyyy 2  -1-  czzz 2 
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Constructing  this  solution  explicitly,  interface  matching  can  then  be  implemented  not  just  vertically  but 
also  horizontally,  avoiding  the  construction  of  artificial  adiabatic  barriers  between  adjoining  sections  of 
surface  metal. 

For  simple  cuboids  with  uniform  surface  fluxes  P1,F2)--)^6  the  solution  reduces  to  just  x,y,z  and 
x2,  y2,  z2  terms,  which  gives  a  simple  formulation  of  a  finite  element  approach  based  on  interface  matching 
of  primitive  cuboids.  This  allows  ready  treatment  of  MMICs  and  of  hybrid  MICs  with  arbitrarily  detailed 
surface  metallisation. 


4.5  MMIC  with  vias 

In  Section  4.1,  a  simple  analytical  expression  was  derived  for  the  thermal  resistance  matrix  of  a  homo¬ 
geneous  MMIC.  In  this  section  it  is  shown  how  to  derive  the  thermal  resistance  matrix  for  a  MMIC 
containing  an  arbitrary  distribution  of  full  or  partial  thickness  vias,  or  with  partial  substrate  thinning. 

Analytical  solution  for  a  single  via  has  been  presented  previously  [58],  but  this  solution  assumed  a 
periodic  distributions  of  vias,  so  that  simple  adiabatic  boundary  conditions  could  be  imposed  on  a  surface 
enclosing  each  via.  For  the  MMIC  and  hybrid  MIC  case  considered  here,  no  such  periodicity  will  generally 
exist. 
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Bonani  et  al  have  presented  approximations  both  for  vias  and  for  partial  substrate  thinning.  Their 
description  of  vias  is  particularly  simple  and  effective,  and  represents  the  MMIC  with  vias  by  means 
of  the  thermal  resistance  matrix  for  the  corresponding  homogeneous  MMIC,  with  parallel  resistances 
constructed  for  stand-alone  metal  vias,  to  take  account  of  the  greatly  increased  thermal  conductivity  of 
the  via  metal  compared  to  the  GaAs  die  [1]— [3]. 

Though  vias  have  been  shown  to  be  largely  unimportant  in  thermal  simulations  [4] ,  (but  not  totally 
insignificant  [1]),  it  may  be  the  case  that  the  thermal  effects  of  vias  could  become  particularly  relevant 
with  joint  developments  in  via  technology  and  development  of  accurate  design  tools.  Partial  substrate 
thinning,  by  contrast,  is  known  to  have  a  major  impact  on  calculated  MMIC  temperatures  [l]-[3]. 

Increases  in  available  computer  memory  and  processor  speed  mean  that  it  is  now  possible  to  construct 
and  utilise  analytical  solutions  which  would  previously  have  been  of  no  practical  value.  The  analytical 
solutions  presented  here  require  inversion  and  diagonalisation  of  large  matrices  (of  order  up  to  several  x  103 
squared).  Such  solutions  are  therefore  vastly  more  computationally  expensive  than  the  simple  resistance 
matrix  approximations  of  Bonani  et  al,  but  they  represent  an  analytical  solution  that  is  numerically  exact 
when  fully  converged,  and  that  becomes  increasingly  more  tractable  with  increased  computing  power. 
The  computational  time  required  to  construct  such  accurate  thermal  resistance  matrix  descriptions  has 
no  impact  on  coupled  electro-thermal  run  times,  as  the  resistance  matrices  are  calculated  prior  to  coupled 
electro-thermal  solution.  The  analytical  solution  for  full  thickness  vias  generalises  to  provide  a  solution  for 
partial  thickness  vias  or  partial  substrate  thinning.  These  solutions,  and  construction  of  the  corresponding 
thermal  resistance  matrices,  are  now  described. 

To  the  authors’  best  knowledge,  this  section  presents  the  first  analytical  solution  for  the  3-dimensional 
temperature  distribution  in  a  cuboid  MMIC  with  an  arbitrary  distribution  of  vias  of  arbitrary  cross-section 
(full  or  partial  thickness). 


4.5.1  Formulation  of  the  linearised  problem 


By  expressing  the  total  thermal  conductivity  as  the  sum  of  GaAs  MMIC  conductivity,  plus  an  additional 
contribution  in  the  vicinity  of  each  via,  j  =  1,...,  A,  the  steady-state  heat  diffusion  equation  for  a  MMIC 
containing  N  vias  can  be  written 


V. 


k(T)  +  ^Hj(x,y)  [A kj(T)  +  8k  j (T)\  l  VT 


(53) 
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I 


where  T  is  physical  temperature;  k(T)  is  the  temperature  dependent  thermal  conductivity  of  the  GaAs 
MMIC;  Akj(T)  for  j  =  are  the  (large)  differences  between  the  thermal  conductivities  of  the 

metal  vias  and  the  GaAs  MMIC,  assuming  the  same  functional  form  for  temperature  dependence  as 
k(T);  Skj(T)  for  j  =  1, N  are  the  small  perturbations  on  A kj(T)  due  to  the  difference  in  temperature 
dependence  of  the  via  metal  thermal  conductivity  from  /c(T);  and 

1  in  the  region  of  via  j 

Hj(x,y)  =  <  .  (54) 

0  otherwise 

Then  neglecting  the  small  quantities  Skj(T ),  this  equation  can  be  written 

V.  [k(T)VT]  = -V.  ^2Hj(x,y)AKj(T)VT  .  (55) 

_  j 

So  the  heat  diffusion  equation  for  the  MMIC  with  vias  has  been  converted  into  an  equation  with  uniform 
GaAs  MMIC  thermal  conductivity  k(T),  but  with  a  temperature  dependent  volume  heat  sink  term, 
described  by  the  right-hand  side  of  Eq.  (55). 

The  justification  for  dropping  the  terms  SKj(T)  is  discussed  in  subsection  4.5.5  below. 

Eq.  (55)  is  non  linear  in  T .  To  linearise  the  equation,  the  Kirchhoff  transformation  is  performed  [55] 
giving 

V20  =  -V. 

with  the  £j  defined  in  subsection  4.5.5  below. 

In  terms  of  linearised  temperature,  6 ,  the  boundary  conditions  become 


where 


dxlx=0'L 

do . 

6{z  =  D) 


0, 
0, 
Ts , 


$0  j  _  Si{x>y)  D 

KSTzU=0  ~  +  ZjtiHj(x,y)  i 


Si{x,y) 


1  for  (x,  y)  €  A 
0  otherwise 


(57) 

(58) 

(59) 

(60) 

(61) 


29 


4.5.2  Double  cosine  transform 


To  solve  this  problem  the  double  Fourier  series  solution  of  Eq.  (14)  is  generalised  to 


6  =  COS  XmX  COS  UnyZmn  0), 


_  TO7T 

Xm  ~  T’ 


and  m,  n  =  0, 1, 2, ... . 

It  is  apparent  that  Eq.  (62)  immediately  satisfies  the  adiabatic  x-  and  y- flux  boundary  conditions 
given  by  Eqns.  (57)  and  (58).  Substituting  the  z-direction  boundary  conditions,  Eqns.  (59)  and  (60), 
into  solution  Eq.  (62)  provides  conditions  on  the 


Z>mn  |  z—D  —  Ts  $0 m  &0 


KsLW  (1  +  <5mo)(l  +  <^no) 


Eqns.  (66)  and  (23)  assume  that  no  elementary  heating  element,  Du  ever  lies  totally  or  partially  over 
a  via.  For  semiconductor  active  device  channels  and  full  thickness  metal  vias  this  will  always  be  true, 
but  the  extension  to  more  general  surface  heating  elements,  partial  thickness  vias  and  partial  substrate 
thinning  is  immediate. 

Finally,  performing  a  double  cosine  transformation  of  Eq.  (56)  [9]  gives 


0  —  ^  y  ^(m'n')(jnn)* 
(mn) 


^(m'n')()nn)  —  ^  (l  "f*  $m/o)  (1  H"  ^n'o)  ^mm'  ^nn' 

py, 

'  /  j  ^J  mm'  nn' 


(m'n')(mn) 


(l  +  ^m'o)  (1  $n'o)  (*  m'  +  Pn'1) 
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+]C& 


(>Zn+ti)Cn'*&> 


(cos  Am/  x f  sin  Am  -  cos  A mt  xf  sin  Xmx j* ) 

+  C^X 

(cos  Ain'yf  Sin/zny^  -  cos  sin  Mn2/f) 


and  where,  for  simplicity,  the  vias  have  been  assumed  to  have  uniform  rectangular  cross-section  with 
faces  parallel  to  the  faces  of  the  MMIC.  xf,xf  are  the  a: -coordinates  of  the  edges  of  via  j  parallel  to  the 
y— axis,  and  yf  are  the  y— coordinates  of  the  edges  of  via  j  parallel  to  the  x— axis.  The  generalisation 
to  vias  of  arbitrary  uniform  cross-section  is  immediate. 

Elementary  integrals  I^m, ,  l^yn,  have  been  defined  as 


rxj 

mm'  =  /  C0sAmXC0sA  m'Xdx, 

Jxf 

fVj 

Inyn>  =  /  cos  n„y  cos  /j,n'y  dy 

Jy? 


and  are  simple  to  evaluate  analytically. 

The  linearised  partial  differential  equation  for  the  GaAs  MMIC  with  vias,  Eq.  (56),  with  boundary 
conditions  Eqns.  (57)-(60),  has  thus  been  converted  to  a  system  of  coupled,  linear,  2nd  order,  ordinary 
differential  equations,  Eq.  (67),  for  the  Zmn  of  Eq.  (62),  with  boundary  conditions  given  by  Eqns.  (65) 
and  (66). 

For  the  special  case  of  —  0  for  all  j,  corresponding  to  a  MMIC  with  no  vias,  this  system  of  equations 


just  reduces  to 


d?  Zmn 
dz2 


■  o'2  Z 

hnn  ^Tnn  y 


where  7 =  A^  +  /i2 ,  with  general  solution 


Zmn  — 


A  +  Bz,  for  (m,  n)  =  (0, 0) 

Cmn  cosh7mnz  -f  Smn  sinh7mn2,  otherwise 


corresponding  to  the  usual  separation  of  variables  solution. 
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4.5.3  Solution  of  the  coupled  linear  system 

In  matrix  notation,  the  coupled  system  of  equations,  Eq.  (67)  can  be  expressed 

^-4-^  =  0-  (74) 

Examination  of  Eqns.  (68)  and  (69)  shows  that,  in  general,  all  elements  of  matrix  ^  will  be  non  zero. 
For  the  case  of  a  MMIC  with  no  vias,  A  will  be  diagonal  with  all  diagonal  elements  non  zero.  However, 
although  the  matrix  elements  of  matrix  £  are  mostly  non  zero  in  the  general  case,  putting  (m,  n)  =  (0, 0) 
in  Eq.  (69)  means  that  #(m'n')( oo)  is  seen  to  be  zero  for  all  (m^n').  For  the  case  of  a  MMIC  with  no 
vias,  B  is  diagonal  with  one  diagonal  element  equal  to  zero. 

Writing  A”1,  J|  and  Z  in  partitioned  form,  with  (double)  index  (m',  n')  =  (0, 0)  corresponding  to  the 
first  row  and  (m,  n)  =  (0, 0)  corresponding  to  the  first  column, 


.4"1  = 


a-i  A-11 
A11  - - 12 


0  Bj2 


e  £» 


Then  partitioning  A  1  B_  in  the  same  way  as  above 


d2Z00/dz2 
d2Z!/dz 2 


0  A~1B\ 


0  A~2Br 


This  gives  rise  to  the  closed  set  of  equations  for  Z_\  he.  for  all  the  Zmn  except  Zq, 

Z'  =  0 

dz2  ^=22-  2, 


subject  to  the  boundary  conditions 


Zl\z-d  —  Q 


and  Eq.  (66)  for  all  m,  n  excluding  (m,n)  =  (0,0). 

Solving  the  homogeneous  linear  system  for  then  gives  the  inhomogeneous  equation  for  Zoo, 


d2Z00  _  !  t  , 

^dz2  ”  —12—  i 


32 


subject  to  the  boundary  conditions 


Zqq  \z=d 


dZoo 

dz 


Ts , 


(82) 

(83) 


The  coupled  system  for  Z7,  Eq.  (79),  can  be  solved  by  following  the  standard  procedure  for  systems 
of  second  order  linear  differential  equations  with  constant  coefficients  [59].  Define 


then  it  is  immediately  obvious  that, 


(84) 


(85) 


where 


Writing 


G  =  e^u, 


(86) 


(87) 


where  from  the  definition  of  Eq.  (84) 


and  u1  is  independent  of  z,  Eq.  (85)  gives  the  eigenvalue  problem 

(4±i22  -72iy =e- 


(88) 


(89) 


For  the  zero  via  case,  the  72  are  pure  real  and  positive  definite,  going  as  A^  +  /z2  with  (m,  n)  ^  (0, 0). 
However,  the  eigenproblem  of  Eq.  (89)  is  constructed  for  a  non  symmetric  real  matrix,  A^1B  22,  so  the 
eigenvalues  can  be  complex.  Although  not  obvious  for  the  most  general  case,  consideration  of  the  terms 
of  the  matrices  A  and  £  of  Eqns.  (68)  and  (69)  shows  that  typically  the  additional  terms  due  to  vias  will 
be  a  small  perturbation  on  the  positive  definite  diagonal  terms  of  A~lB ^  corresponding  to  A^  +  /i2 .  It 
is  therefore  to  be  expected  that  the  modified  eigenvalues  of  A~lB will  remain  real  and  positive  in  the 
presence  of  vias.  This  is  found  to  be  the  case  for  all  calculations  performed  so  far.  In  the  general  case 
det(Al^g22)  will  be  non  zero,  so  that  7  is  never  zero  (for  non  trivial  solutions  u1). 
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Solving  Eq.  (89)  by  standard  numerical  techniques,  to  obtain  eigenvalues  7?  and  eigenvectors  uj,  the 
general  solution  to  Eq.  (85)  can  be  written 

z'  =  J2  (ste+7i* + s'e”7jZ)  4*  (90) 

3 

To  obtain  the  unknown  coefficients  cf,  the  boundary  conditions  given  by  Eqns.  (66)  and  (80)  are  imposed. 
The  vectors  c*  of  components  cf  are  then  seen  to  satisfy 

£^(c+-<r)  =  g>  (91) 

C/+c+  +  C/_c”  =  0,  (92) 


where  the  columns,  of  matrix  J77  are  given  by 


similarly 


Kj.,  =  e+7j£V,-, 


and  p  is  the  vector  with  components 


e-7'V, 


P(mn)  “  ksLW(1  +  Smo)(l+5n0)  Y^P™Pi 

for  (m,n)  7^  (0,0). 


Having  constructed  the  solution  for  Z'  the  solution  for  Z0o  is  immediately  obtained  as 


Zoo  =  A  +  Bz  +  Y, 4  (S+e+7j2  +  cj e-7i*) , 


and  arbitrary  constants  are  determined  by  imposing  boundary  conditions,  Eqns.  (82)  and  (83). 

This  gives  a  complete  analytical  solution  for  the  case  of  a  cuboid  MMIC  with  an  arbitrary  distribution 
of  vias.  The  major  computational  limitation  on  this  approach  is  the  order  of  the  matrix  problems  required 
for  full  convergence  of  the  double  Fourier  series  expansion,  Eq.  (62).  For  a  MMIC  with  no  vias,  adequate 
convergence  is  typically  found  to  be  obtained  for  m,  n  =  0, ~  50,  dependent  on  the  size  of  the  smallest 
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feature  to  be  resolved  relative  to  the  die  size  [10,  11].  The  inclusion  of  vias  is  often  found  to  be  a  small 
perturbation  on  the  MMIC  solution  [1,  4],  so  similar  convergence  rates  might  be  expected  for  the  solution 
with  vias.  However,  the  effects  of  discontinuities  in  thermal  conductivities  is  to  introduce  discontinuities 
into  thermal  gradients  at  the  via  interfaces.  These  gradient  discontinuities  could  require  larger  number 
of  terms  for  an  accurately  converged  solution. 

Assuming  adequate  convergence  of  Eq.  (62)  with  m2  +  n2  <  502  for  the  system  under  consideration, 
the  size  of  the  matrices  to  be  constructed  is  of  order  ~  2000  x  2000,  requiring  ~32  MB  of  real  double 
precision  storage  for  each  matrix.  Inversion,  diagonalisation  and  multiplication  of  matrices  of  this  size 
typically  takes  ~1  hour  to  ~10  hours. 


4.5.4  Thermal  resistance  matrix 


To  construct  the  thermal  resistance  matrix,  analytical  expressions  are  required  for  the  average  tempera¬ 
tures  over  elementary  heating  elements  £>*.  From  Eqns.  (62)  and  (23)  these  are  given  by 


SJm^\z=odxdy 

IJDidxdy 


(99) 

(100) 


or  in  matrix  notation, 


eiv  =  zQO\o+LiT  z!\0, 


(101) 


where 


fi  _  1mn 

J(mn)  ~  ji  > 
i00 


(m,n)  ^  (0,0). 


Prom  Eqns.  (90),  (91)  and  (92), 


l'|o  =  flp, 


where 

=  (l+lLz.~1u±) 

and  where  U_  is  the  matrix  with  columns  U_j  =  .  Hence, 

&av  =  ^oolo+f  TS2 


(102) 


(103) 


(104) 


(105) 


and 


@av  ~  Z00\o  =  ^2/RTHikPk, 
k 


(106) 
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where  the  thermal  resistance  matrix  __  is  given  by 

—1  £1 

RTHik=LiTRL,  (107) 

with 

~  k  —4  k 

I(mn)  =  KS LW  (l  +  5m0)(l+<5no)/mn'  (108) 

As  for  the  homogeneous  case,  the  thermal  resistance  matrix  is  seen  to  be  obtained  purely  in  terms 

of  material  parameters,  and  the  geometry  of  the  MMIC.  As  it  is  derived  for  the  linearised  temperature, 
it  is  independent  of  both  temperature  and  power  density,  and  hence  of  electrical  bias  point  (though 
parameters  must  be  chosen  to  represent  typical  operating  conditions).  The  matrix  only  has  to  be 
constructed  once,  for  repeated  use  in  coupled  electro-thermal  simulations,  and  whatever  the  complexity 
of  the  MMIC  internal  structure  the  order  of  the  thermal  resistance  matrix  is  determined  purely  by  the 
(small)  number  of  surface  heating  elements,  A- 

4.5.5  Non  linear  thermal  conductivity  perturbation 

In  subsection  4.5.1,  terms  Skj(T)  were  dropped  from  the  heat  diffusion  equation  on  the  grounds  that 
they  were  a  small  perturbation  on  terms  A Kj(T).  This  assertion  is  now  justified.  Also,  in  the  linearised 
Eq.  (56),  the  constants  £j  were  introduced;  these  are  now  defined. 

If  the  thermal  conductivity  in  the  GaAs  MMIC  body,  region  I,  is  k/(T),  and  the  thermal  conductivity 
in  the  metal  vias,  region  II,  is  K//(T),  (and  for  simplicity  all  vias  are  assumed  to  be  made  of  identical 
material,  though  this  is  not  a  necessary  restriction),  total  thermal  conductivity  throughout  the  whole 
MMIC,  Ki0t(T),  can  be  written 

ntot(T)  =  kj(T)  +  J2Hj(x,y)  [kji(T)  -  k/(T)]  .  (109) 

3 

To  a  good  approximation,  the  thermal  conductivity  of  GaAs  can  be  assumed  to  have  a  simple  power  law 
dependence  on  physical  temperature  over  temperature  ranges  of  interest  for  device  operation  [56],  and 
the  thermal  conductivity  of  via  metal  can  be  assumed  independent  of  temperature  [T] , 


KlCO  = 

(T\ 
KI  \TsJ 

-fa 

|  b  =  1.22  for  GaAs, 

(110) 

kii{T)  = 

K II 

for  via  metal. 

(111) 

Then  the  A k,-(T)  and  8Kj(T)  for  each  via. 

,  3  =  !»•• 

TV,  can  be  defined  by  the  relations, 

nn(T)  -  *j(T)  =  A k5{T)  Hh  ^(T),  (112) 


36 


A kj(T)  ==  (rjjKii  -  «/) 


—  £jKi(T)i 


which  defines  the  £j,  and 


Shij(T)  =  ku  -  (1  +  £j)  ki(T). 


A Kj(T)  of  Eq.  (113)  equals  the  difference,  ku{T)  -  kj(T),  assuming  the  same  functional  form  for  the 
temperature  dependences  of  k n(T)  and  kj(T).  Eq.  (114)  then  gives  the  small  perturbation  on  A Kj(T) 
resulting  from  the  fact  that  the  temperature  dependences  of  kjj(T)  and  Ki(T)  are  actually  of  different 
functional  form.  This  separation  of  A/c j(T)  and  Skj(T)  contributions  to  kh(T)  -  «/(T)  is  necessary  to 
allow  linearisation  of  Eq.  (55)  by  use  of  the  Kirchhoff  transformation,  Eq.  (9). 

The  constants  rjj  (and  hence  £j)  in  Eq.  (113)  are  introduced  to  allow  minimisation  of  the  errors 
introduced  by  dropping  terms  6kj(T)  from  Eq.  (53).  Although  ku{T)  -  «j(T)  is  the  same  in  each  via, 
these  constants  are  given  different  values  in  each  via,  i ]j  (and  so  £j)  j  =  l,...,iV,  to  allow  for  the  fact 
that  the  temperature  distribution  will  differ  from  via  to  via.  The  use  of  the  £j ,  in  minimisation  of  the 
error  due  to  neglect  of  differences  between  the  functional  forms  of  the  temperature  dependences  of  the 
thermal  conductivities,  is  now  discussed.  The  same  general  approach  applies  to  hierarchical  construction 
of  thermal  resistance  matrices  for  metal  on  MMIC  structures,  as  described  in  Section  4.4. 

The  magnitude  of  the  perturbation  Shij(T)  on  A kj(T)  is  given  by 

Ski(T)  =  1  \lLL  ( IS)  _(1  +  £.)1  (115) 

A «j(T)  &  [«/  \TsJ  [ 

This  is  seen  to  be  exactly  zero  for  the  case  that  the  square  bracket  on  the  right  hand  side  of  Eq.  (115) 


goes  to  zero.  So  defining  the  £j  by 


i  +  ij  — 


Ku  (Tx 


V  Ts 


for  typical  vias  temperatures  Tvia. ,  and  writing  the  temperature  variation  over  via  j  a =  Tviaj  ±  A  Tj} 

ism,, _ ! _ di7) 

Attj(T)  _ kj  (  ts  \  \  Tviaj ) 

KII  J 


which  is  given  approximately  by 


8  k  AT) 

A  k^T) 


A  Tj 
t  ±1.4— — — 


for  a  typical  ratio  kj/ku  of  ~  1/8. 
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For  a  variation  across  a  via  of  ±10  K,  from  a  typical  via  temperature  of  ~300  K,  this  error  is 
6Kj(T)/AKj(T)  «  ±5%.  The  impact  on  calculated  active  device  temperatures,  due  to  neglect  of  errors 
of  this  magnitude  in  the  via  thermal  conductivities,  is  expected  to  be  small. 

Clearly,  mean  via  temperatures  Tvia.  are  not  known  a  priori ,  so  they  have  to  be  obtained  self- 
consistently,  for  instance  by  relaxation  starting  with  temperature  values  obtained  from  the  simple  ana¬ 
lytical  solution  for  the  MMIC  without  vias.  This  self-consistent  solution  for  the  £j  then  gives  an  approxi¬ 
mate  treatment  of  the  non  linear  perturbation  arising  from  the  difference  in  functional  form  between  the 
temperature  dependences  of  the  thermal  conductivities  in  the  GaAs  MMIC  and  in  the  metal  vias. 

As  the  metal  vias  are  of  small  cross-sectional  area,  and  have  high  thermal  conductivities,  the  greatest 
variation  in  temperature  is  expected  to  be  along  the  2-axis,  corresponding  to  the  heat  flux  from  the  active 
devices  at  the  MMIC  surface  to  the  heat  sink  mount  at  the  MMIC  base.  The  accuracy  of  the  solution  can 
therefore  be  improved  further,  by  subdividing  the  vias  in  the  2-direction.  By  generalising  the  mean  via 
temperature  description,  Tviaj,  and  solving  self-consistent ly  for  mean  via  temperatures,  Tviajk ,  in  each 
vertical  element,  fc,  of  any  given  via,  j,  typical  temperature  variations  across  via  elements,  A Tjk>  can  be 
reduced.  The  description  of  the  vias,  j  =  1, IV,  already  allows  for  the  possibility  that  the  cross-section 
of  each  physical  via  can  be  constructed  from  mathematical  sub- vias  of  smaller  cross-section.  Combined 
with  the  generalisation  to  vertical  subelements,  fc,  described  by  mean  temperatures  TViajky  this  means 
that  the  variation  of  temperatures  across  elementary  via  volumes,  A Tjk>  can  be  made  arbitrarily  small 
giving  an  essentially  numerically  exact  self-consistent  solution  to  the  full  non  linear  problem  of  a  MMIC 
with  vias.  However,  this  full  solution  is  very  computationally  demanding  and  generally  not  warranted  by 
the  degree  of  inaccuracy  due  to  the  non  linear  perturbation. 

In  a  fully  coupled  electro-thermal  solution,  it  would  be  impractical  to  obtain  the  fully  non  linear 
solution  for  the  MMIC  thermal  resistance  matrix  at  each  iteration  of  a  global  electro-thermal  solution. 
However,  by  obtaining  the  thermal  resistance  matrix  for  a  typical  operating  point,  the  resulting  linear 
description  of  the  MMIC  should  be  sufficiently  accurate  for  description  of  the  global  system  over  a 
satisfactory  range  of  operation. 

Figure  9  shows  the  calculated  temperature  distribution  of  a  GaAs  MMIC  bearing  two,  6-finger,  power 
transistors.  Figure  10  shows  the  same  die  with  a  central  via.  For  speed  of  calculation,  these  illustrative 
simulations  used  just  a  small  number  of  basis  states  given  by  the  relation  m2  +  n2  <  202.  The  details 
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of  the  individual  device  fingers  are  therefore  not  resolved,  but  the  impact  of  the  centrally  placed  via  is 
clear,  reducing  central  temperature  and  modifying  the  temperatures  in  the  vicinity  of  the  active  device 
fingers. 


4.5.6  Partial  thickness  vias  and  substrate  thinning 


To  treat  the  case  of  partial  thickness  vias  and  partial  substrate  thinning,  and  to  allow  more  accurate 
minimisation  of  non  linear  perturbations,  6nj(T),  the  extension  to  vias  discretised  along  the  z-axis  is  now 
presented.  The  linearised  heat  diffusion  equation,  Eq.  (56),  becomes 


V20  =  -V. 


'£tikHJ{x,y)Hk(z)Ve  , 

jk 


(119) 


where 

( 

1,  within  MMIC  layer  k 
0,  otherwise 

Then  the  coupled  system  for  the  Zmn  of  Eq.  (67)  becomes 


Hk(z)  =  { 


(120) 


d_ 

dz 


-R(z)Z  =  0, 


(121) 


where  A  and  B  are  now  piecewise  constant  and  have  the  same  form  as  previously  within  each  MMIC 
layer,  Eqns.  (68)  and  (69),  but  with  -4  Yljk  ZjkHk{z).  Integrating  Eq.  (121)  across  an  interface, 

it  is  seen  to  reduce  to  the  original  Eq.  (67)  within  each  layer,  with  continuity  of  J[  and  A dZ_/dz  across 
each  interface.  This  problem  can  then  be  treated  by  a  transfer  matrix  approach. 

The  above  formulation  can  also,  in  principle,  treat  the  case  of  totally  arbitrary  surface  metallisation 
considered  as  an  effective  partial  thickness  via  of  complicated  cross-section.  However,  treatment  of  surface 
detail  that  is  too  fine  will  generally  produce  convergence  difficulties.  Also,  this  approach  is  more  limited 
in  its  description  of  non  linear  interface  matching  between  metal  and  GaAs. 


4.5.7  Alternative  formulation 

After  some  manipulation,  Eq.  (56)  can  be  written  in  the  equivalent  form, 

V20  =  X>(1  +  Zi)  VHj(a:,y).Vfl,  (122) 

3 

as  long  as  £?  >  -1  for  all  j.  This  form  again  gives  rise  to  Eq.  (67)  with  matrix  £  of  similar  form  to  Eq. 
(69),  but  now  matrix  A  is  diagonal  and  is  trivial  to  invert.  This  represents  a  significant  saving  in  the 
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construction  of  A  1JB  when  the  order  of  matrix  A  is  large.  A  is  given  by 
A-1  jE?(m/nq(mn)  =  (A^/  +  Mn')  drum' dun'  + 

? ln  (1  +  ^  LW(l  +  <Wo)(l  +  <5n'0) 

3 

K,  +  «0  /£».«. 

+  Inn’^m  X 

(cos  Am/a;j'sinAm:rj'  -  cos  Am/ xj1  sin  \mxf) 
x 

+  Imm'tJ'nX 

(cos  Hn’yf  sin  jinyf  -cosn„>yf  sin/vyj7) 

(123) 

4.5.8  Approximate  treatment  of  vias 

The  double  Fourier  series  solutions  for  the  treatment  of  vias,  presented  above,  are  computationally 
demanding.  However,  the  simple  equivalence  principle  approach  of  Bonani  et  al  [1]— [3]  is  conceptually 
simple  and  is  cheap  to  implement.  It  is  described  here  within  the  framework  of  the  resistance  matrix 
formulation. 

Partitioning  the  thermal  resistance  matrix  equation,  Eq.  (7),  by  elementary  surface  elements  corre¬ 
sponding  to  active  devices,  A 0a,£_a  and  via  holes,  A#V,.F  gives 


with 

M«  =  Rth(-Ev)  (125) 

describing  the  fluxes  through  the  parallel  equivalent  via  thermal  resistances,  Rth- 

All  vias  are  assumed  identical  for  simplicity.  The  negative  sign  in  Eq.  (125)  reflects  the  fact  that 
the  parallel  resistances,  corresponding  to  the  greater  thermal  conductivity  of  the  metal  compared  to 
GaAs  and  applying  the  equivalence  principle,  are  all  effectively  external  to  the  GaAs  die.  (This  means 
that  although  the  parallel  resistances  describe  increased  heat  flow  through  via  metal  in  an  approximate 
fashion,  they  do  not  correctly  include  the  effects  of  increased  heat  flow  through  the  side  walls  of  the  via 
metal  in  the  GaAs.) 
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Combining  Eqns.  (124)  and  (125),  the  global  thermal  resistance  matrix  relating  active  device  tem¬ 


perature  rises  to  active  device  power  dissipations  is  given  by, 


Ma  =  R9TH^ 


(126) 


with 

Eth  =  Eth  -  Eth  ( Eth  +  L RTH )  _1  RVth-  (!27) 

For  a  via  of  uniform  cross-section,  passing  through  a  GaAs  die  of  thickness,  D,  with  difference  between 
metal  and  GaAs  thermal  conductivity  given  by,  Ak,  the  via  thermal  resistance  is  given  simply  by 


Tzvia  _  ^ 


(128) 


To  describe  more  accurately  the  thermal  non  linearity  due  to  the  difference  in  functional  form  between 
the  temperature  dependent  thermal  conductivities  of  metal  and  GaAs,  An  can  be  ascribed  a  different 
value  in  each  via,  corresponding  to  respective  mean  temperature. 

Combining  this  approach  with  a  simple  generalisation  of  the  2-layer  thermal  resistance  matrix  solu¬ 
tion  of  Section  4.3  provides  a  straightforward  description  of  partial  thickness  vias  and  partial  substrate 
thinning.  In  the  case  of  partial  substrate  thinning  or  ‘bathtubs’  [1],  flux  at  the  surface  of  the  partial 
thickness  bathtub  metallisation  is  assumed  piecewise  constant,  rather  than  completely  uniform. 


4.6  Large  area  substrates 


For  large  area  substrates,  from  MMIC  grid  arrays  up  to  circuit  board  level,  surface  fluxes  due  to  convection 
and  radiation  must  be  included  in  the  solution  of  the  heat  diffusion  equation.  (The  same  would  hold  for 
small  area  substrates  under  any  circumstances  in  which  convection  became  important.)  For  moderately 
low  power  densities  and  temperatures  these  fluxes  can  be  described  by  the  linear  radiation  boundary 
condition,  Eq.  (13)  [19].  Putting  a0,D  =  1,  H0jD  ±  0  the  steady-state  heat  diffusion  equation,  Eq.  (10), 
can  be  solved  to  obtain  the  corresponding  thermal  resistance  matrix  for  a  substrate  with  radiative  and 
convective  surface  losses.  Figures  11-13  show  the  results  of  electro-thermal  simulations  performed  using 
such  a  thermal  resistance  matrix. 

Putting  a0  =  aD  =  1,- H0  -  HD  =  H,0o(x,y)  =  0D(x,y)  =  Oo,po{x,y)  =  0  and  p0(x,y)  = 
£\  Si{x,y)Pi,  the  analytical  solution  Eq.  (14)  is  obtained  with  coefficients 


A  = 


(129) 
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B 


B  = 

(130) 

Cmn  — 

H  tanh-ymnD  +  Ko-fmn 

mn  H  +  K07mn  tanh  7mn-D  ’ 

(131) 

Bmn  “ 

H-j-KQ'yrnn  1 3. fill  'Ymn  D 

2ffKo7m„+(^7jin+H2)  tenh 7 mnD 

(132) 

X  £W(l+«m0)(l+«no)  Ann-P*) 

where  k0  is  the  thermal  conductivity  of  the  substrate. 

The  corresponding  expression  for  the  thermal  resistance  matrix,  describing  the  temperature  rise  of 
the  surface  heating  elements  above  Oq  ,  is  given  by 


RrHi}  = 


(S+2) 


\  p 

HLW  00 


+£ 


_ ijcu.ui  imn ~  n,Q  Jmn 

2Hft07mn  +  («o7mn  +  #2)tanh7i 
4 


X 


LW  (1  +  8mQ  )  ( 1  -f  <^no) 


p  p 

•Lmn-Lrr 


i00 


(133) 


Based  on  these  expressions,  the  simulation  of  Figure  11  shows  the  simulated  reverse  side  of  a  Cu/FR-4 
mount  bearing  a  3-stage  balanced  amplifier  MMIC.  The  simulation  used  a  two-level  structure  with  thermal 
resistance  matrices  constructed  for  the  partial  Cu  layer  and  the  FR-4  mount  and  matched  at  the  interface. 


(This  construction  was  for  illustrative  purposes  only.  An  analytical  solution  can  be  constructed  directly 
for  this  configuration,  consisting  of  a  double  Fourier  series  in  each  layer,  with  Fourier  components  fully 
matched  by  solution  of  a  corresponding  linear  system.)  Two-layer  solutions  of  this  type  require  imposition 
of  the  radiation  boundary  condition  on  the  surfaces  of  both  layers,  so  as  to  produce  exact  cancellation  of 
‘fictitious’  radiative  and  convective  fluxes  at  the  matched  interfaces. 

As  the  electro-thermal  model  currently  contains  no  explicit  description  of  external  fluid  flow,  the 
surface  flux  parameters  JT0,d  in  the  radiation  boundary  condition,  Eq.  (13),  were  used  as  a  fitting 
parameter  to  experimental  data  (setting  -  H0  =  HD  =  H)  [19].  Only  the  average  temperature  on  the 
back  of  the  Cu  layer  was  matched  to  the  substrate. 

Figures  12  and  13  show  the  impact  of  parameter  variation  on  the  calculated  solution.  Removing  the 
partial  Cu  layer,  Figure  12,  leads  to  pronounced  localisation  of  the  temperature  profile.  The  2  fim  thick 
Cu  layer  therefore  acts  as  an  efficient  heat  spreader  on  FR-4.  Figure  13  shows  the  effect  of  altering  the 
value  of  the  FR-4  thermal  conductivity  to  that  of  AIN  (three  orders  of  magnitude  higher).  The  high  AIN 
conductivity  produces  an  almost  uniform  temperature  across  the  surface  of  the  mount. 
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4.6.1  External  fluid  flow 


The  need  to  use  the  heat  transfer  coefficient  H0,d  as  a  free  parameter,  limits  the  physicality  of  the 
description  provided  by  the  coupled  electro-thermal  model  in  the  case  of  large  area  surface  losses.  For 
larger  area  substrates  #0,0  can  be  obtained  using  standard  convection  correlations  from  the  literature 
[11].  However,  for  the  description  to  be  fully  physical  H0,d  must  be  obtained  from  a  physical  model.  The 
full  3-dimensional  coupled  conduction-convection  problem  is  complex,  e.g.  [60],  however  for  simple  cases 
of  interest,  such  as  laminar  flow  over  a  vertical  plate,  analytical  solutions  can  be  obtained  for  temperature 
distributions  along  the  plate  [61].  This  offers  the  possibility  of  fully  physical  thermal  resistance  matrix 
calculations  to  make  preliminary  studies  of  external  convection  on  device  and  circuit  performance. 

4.7  Time  dependent  temperature 

The  thermal  resistance  matrix  approach  described  so  far  applies  only  to  the  temperature  steady-state. 
However,  the  method  generalises  readily  to  treat  the  thermal  time-dependent  case  in  Laplace  transform 
s-space  [19,  24].  Analytically  generated  thermal  impedance  matrices,  describing  thermal  transients  and 
capacitances,  can  be  used  directly  in  s-space,  for  instance  in  harmonic  balance  simulation,  or  in  the 
time-domain  where  they  correspond  to  the  impulse  response  of  the  thermal  system. 

In  coupled  electro-thermal  simulations  thermal  impedance  matrices  have  to  be  calculated  at  a  small 
number  of  points  in  s-space.  A  modest  amount  of  precomputation  then  determines  all  thermal  time 
evolution  information.  The  coupled  electro-thermal  solution  requires  generalisation  of  the  double  Fourier 
series  approach  to  treat  arbitrary  initial  conditions,  allowing  repeated  resetting  of  initial  conditions  in 
a  time-stepping  solution.  (This  generalisation  also  provides  a  double  Fourier  series  treatment  of  volume 
heat  sources  and  sinks,  without  the  use  of  Green’s  function  techniques.)  The  time-dependent  coupled 
electro-thermal  solution  then  reduces  to  a  series  of  small  non  linear  problems  with  order  given  by  a  small 
multiple  of  the  number  of  active  device  elements. 

The  resulting  formulation  provides  a  fully  coupled,  physical  electro-thermal  solution  in  which  the  main 
approximations  are  flux  and  temperature  discretisation  at  interfaces,  interpolation  of  electrical  power 
dissipation  between  discrete  time-points,  and  numerical  inversion  of  Laplace  transforms.  All  surface  and 
material  non  linearities,  including  temperature  dependent  diffusivity,  are  fully  treated.  This  solution 
scheme  compares  favourably  with  time-stepping  finite  difference  or  finite  element  solutions,  capable  of 


43 


comparable  levels  of  detail  and  accuracy,  which  must  iterate  physical  electrical  and  thermal  solutions  at 
each  time  step. 

5  Summary  and  Conclusions 

The  problem  of  physical,  coupled  electro-thermal  modelling  for  electro-thermal  design  of  high  power 
devices  and  circuits  has  been  introduced.  The  combination  of  the  quasi-2D  Leeds  Physical  Model  with 
the  analytical  thermal  resistance  matrix  approach  has  been  proposed  as  a  viable  solution.  The  Leeds 
Physical  Model  has  been  described  and  the  details  of  the  coupled  electro-thermal  problem  outlined. 
The  motivation  for  the  thermal  resistance  matrix  approach  to  the  thermal  steady-state  case  has  been 
presented. 

The  most  general  formulation  of  the  resistance  matrix  approach  offers  a  potentially  competitive  so¬ 
lution  scheme  for  any  physical  system,  described  by  an  equation  in  some  complex  volume,  which  can 
be  resolved  into  subsystems  allowing  linear  analytical  solution  with  specified  boundary  conditions.  The 
resistance  matrix  approach  then  amounts  to  matching  of  field,  ip,  and  its  derivative,  ip',  at  subsystem 
interfaces,  with  the  resistance  matrix  relating  linearly  ip  and  ip'  at  distinct  interfaces  of  any  given  sub¬ 
volume.  This  approach  produces  a  minimal  system  to  be  solved  for  (discretised)  interface  unknowns, 
thus  providing  full  3-dimensional  solutions  for  the  whole  complex  volume.  It  is  of  particular  value  for 
describing  distinct  physical  systems  coupled  in  a  relatively  small  number  of  isolated  regions. 

Applying  this  approach  to  the  coupled  electro-thermal  system,  the  construction  of  thermal  resistance 
matrices  has  been  described  for  the  case  of  the  homogeneous  MMIC,  bearing  an  arbitrary  distribution  of 
power  transistor  fingers.  It  was  shown  that  this  construction  leads  to  a  small  non  linear  problem  for  the 
coupled  electrical  description,  readily  soluble  by  Newton  methods.  Construction  of  resistance  matrices 
was  then  outlined  for  multi-layer  MMICs  and  MMICs  with  superstructure  such  as  surface  metallisation 
and  air  bridges.  Constructions  were  presented  for  MMICs  with  vias  and  with  partial  substrate  thinning, 
and  the  application  of  the  radiation  boundary  condition  for  larger  area  substrates  was  described. 

The  steady-state  thermal  resistance  matrix  approach  described  here,  only  requires  construction  of  the 
global  thermal  resistance  matrix  for  a  complex  3-dimensional  system  once,  prior  to  the  coupled  electro¬ 
thermal  simulation.  Precomputation  to  construct  the  thermal  resistance  matrices  is  modest,  as  matrices 
for  thermal  subsystems  are  given  by  simple  analytical  expressions.  The  global  thermal  resistance  matrix 
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is  obtained  by  simple  matrix  manipulations  on  the  subsystem  matrices,  to  impose  temperature  and  flux 
matching.  Construction  time  for  the  thermal  resistance  matrices  has  no  impact  on  coupled  electro-thermal 
run  time.  The  magnitude  of  the  non  linear  coupled  electrical  problems  to  be  solved  self-consistently  is 
typically  given  by  the  number  of  active  device  elements,  not  determined  by  the  internal  complexity  of 
the  device  structure. 

Full  treatment  of  thermal  non  linearities  has  been  described,  including  treatment  of  temperature 
dependent  conductivity  (perhaps  requiring  non  linear  interface  matching)  and  non  linear  surface  fluxes. 
Full  treatment  of  interface  matching  and  surface  flux  non  linearities,  based  on  the  linear  thermal  resistance 
matrix  approach,  alters  the  size  of  the  corresponding  coupled  problem,  but  the  size  of  the  non  linear 
system  to  be  solved  can  be  minimised  by  careful  choice  of  surface  elementary  areas,  to  take  full  account 
of  expected  temperature  distributions  and  minimise  the  degree  of  surface  discretisation. 

The  combination  of  the  Leeds  Physical  Model  and  the  analytical  thermal  resistance  matrix  approach 
offers  order  of  magnitude  speed  up  on  coupled  electro-thermal  solutions  based  on  intensive  numerical 
electrical  and  thermal  solutions,  capable  of  comparable  levels  of  detail  and  accuracy.  The  fully  physical 
model  will  provide  a  predictive  tool  for  electro-thermal  design  studies  on  CAD  timescales. 
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7  Captions 


Figure  1:  Self-consistent  Schrodinger-Poisson  solution  from  the  LPM  for  HEMT  band-edge  profile,  sub¬ 
band  structure,  and  carrier  density  normal  to  the  heterointerface  [51]. 

Figure  2:  Typical  FET  cross-section  and  layer  structure  illustrating  Gaussian  boxes  used  in  the  LPM  to 
model  the  2-dimensional  profile  [22]. 

Figure  3:  LPM  time  and  frequency  domain  response  to  a  sinusoidal  5  GHz,  0.45  V  input  signal  [53]. 
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Figure  4:  Thermal  equivalent  circuit  corresponding  to  resistance  matrix  Rij  [18]. 


Figure  5:  Illustrative  GaAs  MMIC  for  generation  of  Rth- 

Figure  6:  Coupled  electro-thermal  calculation  of  surface  temperature  distribution  for  a  60-finger  power 
HEMT. 


Figure  7:  Illustration  of  5x5  power  FET  array  simulated  electro-thermally. 

Figure  8:  Electro-thermal  simulation  of  the  central  power  FET  of  a  5x5  power  FET  array.  Non  linear 
temperature  matching  was  used  at  MMIC-substrate  interfaces. 

Figure  9:  Simulation  of  a  500  pm  x  250  pm  GaAs  MMIC  bearing  two  6-finger  power  transistors  dissipating 
0.3  W  total. 

Figure  10:  Simulation  of  same  MMIC  as  Figure  9  but  with  a  50  pm  x  50  pm  central  via.  15°C  above 
heat  sink  mount  (peak)  to  5°C  above  mount  (edge). 

Figure  11:  Electro-thermal  simulation  of  Cu/FR-4  mount  reverse  bearing  a  balanced  amplifier  MMIC. 
Cooling  entirely  by  radiation  and  convection  (no  heatsink).  46° C  above  ambient  (centre)  to  ambient 
(edge)  [18]  ,[19]. 


Figure  12:  Simulated  temperature  of  the  FR-4  mount  reverse  side  of  Figure  11  neglecting  Cu  layer.  96°C 
above  ambient  (centre)  to  ambient  (edge)  [19]. 

Figure  13:  Simulated  temperature  of  the  FR-4  mount  reverse  side  of  Figure  11  assuming  AIN  parameters. 
10°C  above  ambient  (centre)  to  9°C  above  ambient  (edge). 


54 


Conduction  Band  Edge  (eV) 


o  o  o  o  o 

k>  o  ho  ^  o> 


Carrier  Density  (cm'3) 


1 —  S’O 


Source  Gate  Drain 


Gaussian  Box 


